Changeset 9511


Ignore:
Timestamp:
Jan 27, 2015, 12:46:45 PM (10 years ago)
Author:
davies
Message:

Adding k_nearest_neighbours > 1 to make_nearestNeighbour_quantity_function

Location:
trunk/anuga_core/source/anuga/utilities
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/utilities/plot_utils.py

    r9386 r9511  
    736736    return near_points[keepers], local_coord[keepers]
    737737
    738 ########################
    739 # TRIANGLE AREAS, WATER VOLUME
     738
    740739def triangle_areas(p, subset=None):
    741740    # Compute areas of triangles in p -- assumes p contains vertex information
    … …  
    764763    return area
    765764
    766 ###
    767765
    768766def water_volume(p, p2, per_unit_area=False, subset=None):
    … …  
    779777   
    780778    # This accounts for how volume is measured in ANUGA
    781     # Compute in 2 steps to reduce precision error (important sometimes)
    782     # Is this really needed?
     779    # Compute in 2 steps to reduce precision error from limited SWW precision
     780    # FIXME: Is this really needed?
    783781    for i in range(l):
    784782        #volume[i]=((p2.stage[i,subset]-p2.elev[subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
    … …  
    10501048
    10511049    # Get function to interpolate quantity onto gridXY_array
    1052     gridXY_array=scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()
     1050    gridXY_array=scipy.array([scipy.concatenate(gridX),
     1051                              scipy.concatenate(gridY)]).transpose()
    10531052    gridXY_array=scipy.ascontiguousarray(gridXY_array)
    10541053
    … …  
    10561055    #basic_nearest_neighbour=False
    10571056    if(k_nearest_neighbours==1):
    1058         index_qFun=scipy.interpolate.NearestNDInterpolator(swwXY,scipy.arange(len(swwX),dtype='int64').transpose())
    1059         gridqInd=index_qFun(gridXY_array)
     1057        index_qFun = scipy.interpolate.NearestNDInterpolator(
     1058            swwXY,
     1059            scipy.arange(len(swwX),dtype='int64').transpose())
     1060        gridqInd = index_qFun(gridXY_array)
    10601061        # Function to do the interpolation
    10611062        def myInterpFun(quantity):
    … …  
    10631064    else:
    10641065        # Combined nearest neighbours and inverse-distance interpolation
    1065         index_qFun=scipy.spatial.cKDTree(swwXY)
    1066         NNInfo=index_qFun.query(gridXY_array,k=k_nearest_neighbours)
     1066        index_qFun = scipy.spatial.cKDTree(swwXY)
     1067        NNInfo = index_qFun.query(gridXY_array, k=k_nearest_neighbours)
    10671068        # Weights for interpolation
    1068         nn_wts=1./(NNInfo[0]+1.0e-100)
    1069         nn_inds=NNInfo[1]
     1069        nn_wts = 1./(NNInfo[0]+1.0e-100)
     1070        nn_inds = NNInfo[1]
    10701071        def myInterpFun(quantity):
    1071             denom=0.
    1072             num=0.
     1072            denom = 0.
     1073            num = 0.
    10731074            for i in range(k_nearest_neighbours):
    1074                 denom+=nn_wts[:,i]
    1075                 num+= quantity[nn_inds[:,i]]*nn_wts[:,i]
     1075                denom += nn_wts[:,i]
     1076                num += quantity[nn_inds[:,i]]*nn_wts[:,i]
    10761077            return (num/denom)
    1077 
    1078 
    10791078
    10801079    if(bounding_polygon is not None):

  • trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py

    r9393 r9511  
    1212        domain,
    1313        threshold_distance = 9.0e+100,
    14         background_value = 9.0e+100):
     14        background_value = 9.0e+100,
     15        k_nearest_neighbours = 1,
     16    ):
    1517    """
    1618    Function which makes another function, which can be used in set_quantity
    … …  
    3234            defining the points used to set the new quantity values
    3335        @param domain -- The ANUGA domain
     36        @param k_nearest_neighbors If > 1, then an inverse-distance-weighted average
     37               of the k nearest neighbours is used
    3438        @param threshold_distance -- Points greater than this distance from
    3539            their nearest quantity_xyValue point are set to background_value
    … …  
    3842    OUTPUTS:
    3943        A function f which can be passed to domain.set_quantity('myQuantity', f)
    40     """
     44
     45    """
     46
    4147    import scipy
    4248    import scipy.interpolate
    43 
    44 
    45     if(len(quantity_xyValueIn.shape)>1):
    46         quantity_xyValue=copy.copy(quantity_xyValueIn) # Pointer, no copy
     49    import scipy.spatial
     50
     51    if(len(quantity_xyValueIn.shape) > 1):
     52        quantity_xyValue = quantity_xyValueIn
    4753    else:
    4854        # Treat the single-point case
    49         quantity_xyValue=copy.copy(quantity_xyValueIn.reshape((1,3)))
     55        quantity_xyValue = quantity_xyValueIn.reshape((1,3))
     56
    5057    # Make a function which gives us the ROW-INDEX of the nearest xy point in
    5158    # quantity_xyValue
    52     quantity_xy_interpolator=scipy.interpolate.NearestNDInterpolator(
    53                                 quantity_xyValue[:,0:2],
    54                                 scipy.arange(len(quantity_xyValue[:,2])))
    55 
    56     #
     59    #quantity_xy_interpolator = scipy.interpolate.NearestNDInterpolator(
     60    #    quantity_xyValue[:,0:2],
     61    #    scipy.arange(len(quantity_xyValue[:,2])))
     62
     63    # Make a function which returns k-nearest-neighbour indices + distances
     64    quantity_xy_interpolator = scipy.spatial.cKDTree(quantity_xyValue[:,0:2])
     65
    5766    # Make a function of x,y which we can pass to domain.set_quantity
    5867    def quant_NN_fun(x,y):
    5968        """
    6069        Function to assign quantity from the nearest point in quantity_xyValue,
    61         UNLESS the point is more than 'threshold_distance' away from the nearest point,
    62         in which case the background friction value is used
    63        
     70        UNLESS the point is more than 'threshold_distance' away from the
     71        nearest point, in which case the background value is used
     72
    6473        """
     74
    6575        import scipy
    6676        import scipy.interpolate
     77        import scipy.spatial
     78
    6779        # Since ANUGA stores x,y internally in non-georeferenced coordinates,
    6880        # we adjust them here
    69         xll=domain.geo_reference.xllcorner
    70         yll=domain.geo_reference.yllcorner
    71         z=scipy.zeros(shape=(len(x), 2))
    72         z[:,0]=x+xll
    73         z[:,1]=y+yll
     81        xll = domain.geo_reference.xllcorner
     82        yll = domain.geo_reference.yllcorner
     83        z = scipy.zeros(shape=(len(x), 2))
     84        z[:,0] = x+xll
     85        z[:,1] = y+yll
     86
    7487        # This will hold the quantity values
    75         quantity_output=x*0. +background_value
     88        quantity_output = x*0. + background_value
    7689        # Compute the index of the nearest-neighbour in quantity_xyValue
    77         q_index=quantity_xy_interpolator(z)
     90        neighbour_data = quantity_xy_interpolator.query(z,
     91            k=k_nearest_neighbours)
     92
    7893        # Next find indices with distance < threshold_distance
    79         dist_lt_thresh=( (z[:,0]-quantity_xyValue[q_index,0])**2 + \
    80                          (z[:,1]-quantity_xyValue[q_index,1])**2 < \
    81                         threshold_distance**2)
    82         dist_lt_thresh=dist_lt_thresh.nonzero()[0]
    83         quantity_output[dist_lt_thresh] =\
    84             quantity_xyValue[q_index[dist_lt_thresh],2]
     94        if(k_nearest_neighbours==1):
     95            dist_lt_thresh = neighbour_data[0] < threshold_distance
     96        else:
     97            dist_lt_thresh = neighbour_data[0][:,0] < threshold_distance
     98
     99        dist_lt_thresh = dist_lt_thresh.nonzero()[0]
     100
     101        # Initialise output
     102        quantity_output = x*0 + background_value
     103
     104        # Interpolate
     105        if len(dist_lt_thresh)>0:
     106            numerator = 0
     107            denominator = 0
     108            for i in range(k_nearest_neighbours):
     109                if(k_nearest_neighbours==1):
     110                    distances = neighbour_data[0][dist_lt_thresh]
     111                    indices = neighbour_data[1][dist_lt_thresh]
     112                else:
     113                    distances = neighbour_data[0][dist_lt_thresh,i]
     114                    indices = neighbour_data[1][dist_lt_thresh,i]
     115               
     116                inverse_distance = 1.0/(distances+1.0e-100)
     117                values = quantity_xyValue[indices,2]
     118                numerator += values*inverse_distance
     119                denominator += inverse_distance
     120           
     121            quantity_output[dist_lt_thresh] = numerator/denominator
     122
    85123        return quantity_output
     124
    86125    # Return the quantity function
    87126    return quant_NN_fun

  • trunk/anuga_core/source/anuga/utilities/test/test_quantity_setting_functions.py

    r9437 r9511  
    4343        """
    4444
    45         boundaryPolygon=[ [minX, minY], [minX, minY+100.], [minX+100., minY+100.], [minX+100., minY]]
    46         anuga.create_mesh_from_regions(boundaryPolygon,
    47                                    boundary_tags={'left': [0],
    48                                                 'top': [1],
    49                                                 'right': [2],
    50                                                 'bottom': [3]},
    51                                    maximum_triangle_area = 1.,
    52                                    minimum_triangle_angle = 28.0,
    53                                    filename = 'test_quantity_setting_functions.msh',
    54                                    interior_regions =[ ],
    55                                    verbose=False)
     45        boundaryPolygon=[ [minX, minY], [minX, minY+100.],
     46            [minX+100., minY+100.], [minX+100., minY]]
     47        anuga.create_mesh_from_regions(
     48            boundaryPolygon,
     49            boundary_tags={'left': [0],
     50                         'top': [1],
     51                         'right': [2],
     52                         'bottom': [3]},
     53            maximum_triangle_area = 1.,
     54            minimum_triangle_angle = 28.0,
     55            filename = 'test_quantity_setting_functions.msh',
     56            interior_regions =[ ],
     57            verbose=False)
    5658
    5759        domain=anuga.create_domain_from_file('test_quantity_setting_functions.msh')
    … …  
    103105        inPts=numpy.vstack([xR,yR,zR]).transpose()
    104106
    105         F=qs.make_nearestNeighbour_quantity_function(inPts, domain,
    106                             threshold_distance = 9.0e+100, background_value = 9.0e+100)
     107        F = qs.make_nearestNeighbour_quantity_function(inPts, domain,
     108                threshold_distance = 9.0e+100, background_value = 9.0e+100)
    107109
    108110        # Test that F evaluated in 'ANUGA coordinates' [lower-left = 0,0] is correct
    … …  
    125127        inPts=numpy.vstack([xR,yR,zR]).transpose()
    126128
    127         F=qs.make_nearestNeighbour_quantity_function(inPts, domain,
    128                             threshold_distance = 6., background_value = 9.0e+100)
     129        F = qs.make_nearestNeighbour_quantity_function(inPts, domain,
     130                threshold_distance = 6., background_value = 9.0e+100)
    129131
    130132        # Test that F evaluated in 'ANUGA coordinates' [lower-left = 0,0] is correct
    … …  
    158160
    159161        # Make a polygon-point pair which we use to set elevation in a 'channel'
    160         trenchPoly=[[minX+40., minY], [minX+40., minY+100.], [minX+60., minY+100.], [minX+60., minY]]
     162        trenchPoly = [[minX+40., minY], [minX+40., minY+100.],
     163            [minX+60., minY+100.], [minX+60., minY]]
    161164
    162165        #################################################################
    163166 
    164167        # This example uses a constant, and a raster, to set the quantity           
    165         F=qs.composite_quantity_setting_function([[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],\
    166                                                 domain)
     168        F=qs.composite_quantity_setting_function(
     169            [[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],
     170            domain)
    167171
    168172        # Points where we test the function
    … …  
    176180        # Find the nearest domain point to the second test point
    177181        # This will have been used in constructing the elevation raster
    178         nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
     182        nearest=((domain.centroid_coordinates[:,0]-3.)**2 +
     183                 (domain.centroid_coordinates[:,1]-20.)**2).argmin()
    179184        nearest_x=domain.centroid_coordinates[nearest,0]
    180185        assert(numpy.allclose(fitted[1],-nearest_x/150.))
    … …  
    182187        #################################################################
    183188 
    184         # This example uses a constant, and a raster, to set the quantity, and applies the min/max bound           
    185         F=qs.composite_quantity_setting_function([[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],\
    186                                                 domain,
    187                                                 clip_range = [[-500., 1.0e+100], [-1.0e+100, 1.0e+100]])
     189        # This example uses a constant, and a raster, to set the quantity, and
     190        # applies the min/max bound           
     191        F=qs.composite_quantity_setting_function(
     192            [[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],
     193            domain,
     194            clip_range = [[-500., 1.0e+100], [-1.0e+100, 1.0e+100]])
    188195
    189196        # Points where we test the function
    … …  
    197204        # Find the nearest domain point to the second test point
    198205        # This will have been used in constructing the elevation raster
    199         nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
    200         nearest_x=domain.centroid_coordinates[nearest,0]
     206        nearest = ((domain.centroid_coordinates[:,0]-3.)**2 +
     207                   (domain.centroid_coordinates[:,1]-20.)**2).argmin()
     208        nearest_x = domain.centroid_coordinates[nearest,0]
     209
    201210        assert(numpy.allclose(fitted[1],-nearest_x/150.))
    202211
    … …  
    206215        def f0(x,y):
    207216            return x/10.
    208         F=qs.composite_quantity_setting_function([[trenchPoly, f0], ['Extent', 'PointData_ElevTest.tif']],\
    209                                                 domain)
    210         fitted=F(testPts_X,testPts_Y)
     217        F = qs.composite_quantity_setting_function(
     218            [[trenchPoly, f0], ['Extent', 'PointData_ElevTest.tif']],
     219            domain)
     220        fitted = F(testPts_X,testPts_Y)
    211221        # Now the fitted value in the trench should be determined by f0
    212222        assert(numpy.allclose(fitted[0],50./10.))
    213223        # The second test point should be as before
    214         nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
    215         nearest_x=domain.centroid_coordinates[nearest,0]
     224        nearest = ((domain.centroid_coordinates[:,0]-3.)**2 +
     225                   (domain.centroid_coordinates[:,1]-20.)**2).argmin()
     226        nearest_x = domain.centroid_coordinates[nearest,0]
    216227        assert(numpy.allclose(fitted[1],-nearest_x/150.))
    217228
    … …  
    219230
    220231        # This example uses 'All' as a polygon
    221         F=qs.composite_quantity_setting_function([['All', f0 ], [None, 'PointData_ElevTest.tif']],\
    222                                                 domain)
     232        F = qs.composite_quantity_setting_function(
     233            [['All', f0 ], [None, 'PointData_ElevTest.tif']],
     234            domain)
    223235        fitted=F(testPts_X,testPts_Y)
    224236        # Now the fitted value in the trench should be determined by f0
    … …  
    229241        # This example should fail
    230242        def should_fail():
    231             F=qs.composite_quantity_setting_function([['All', f0], ['All', 'PointData_ElevTest.tif']],\
    232                                                     domain)
     243            F=qs.composite_quantity_setting_function(
     244                [['All', f0], ['All', 'PointData_ElevTest.tif']],
     245                domain)
    233246            # Need to call it to get the error
    234247            F(numpy.array([3.]), numpy.array([3.]))
    … …  
    239252        # This example should fail (since the clip_range minimum >= maximum)
    240253        def should_fail_1():
    241             F=qs.composite_quantity_setting_function([[trenchPoly, f0], ['All', 'PointData_ElevTest.tif']],\
    242                                                     domain,
    243                                                     clip_range = [[ -500., -1000.], [-1.0e+100, 1.0e+100]])
     254            F=qs.composite_quantity_setting_function(
     255                [[trenchPoly, f0], ['All', 'PointData_ElevTest.tif']],
     256                domain,
     257                clip_range = [[ -500., -1000.], [-1.0e+100, 1.0e+100]])
    244258            return
    245259        self.assertRaises(Exception, lambda: should_fail_1())
    … …  
    251265        #
    252266        #
    253         domain=self.create_domain(1.0, 0.0)
     267        domain = self.create_domain(1.0, 0.0)
    254268
    255269        # Evolve the model
    … …  
    259273        # Make a raster from the elevation data
    260274        from anuga.utilities import plot_utils as util
    261         xs=domain.centroid_coordinates[:,0]+domain.geo_reference.xllcorner
    262         ys=domain.centroid_coordinates[:,1]+domain.geo_reference.yllcorner
    263         elev=domain.quantities['elevation'].centroid_values
    264 
    265         allDat=numpy.vstack([xs,ys,elev]).transpose()
    266         util.Make_Geotif(allDat, output_quantities=['ElevTest'], EPSG_CODE=32756,
    267                         output_dir='.', CellSize=1.,k_nearest_neighbours=1)
     275        xs = domain.centroid_coordinates[:,0]+domain.geo_reference.xllcorner
     276        ys = domain.centroid_coordinates[:,1]+domain.geo_reference.yllcorner
     277        elev = domain.quantities['elevation'].centroid_values
     278
     279        allDat = numpy.vstack([xs,ys,elev]).transpose()
     280        util.Make_Geotif(allDat, output_quantities=['ElevTest'],
     281            EPSG_CODE=32756, output_dir='.', CellSize=1.,k_nearest_neighbours=1)
    268282
    269283        # Make a polygon-point pair which we use to set elevation in a 'channel'
    270         trenchPoly=[[minX+40., minY], [minX+40., minY+100.], [minX+60., minY+100.], [minX+60., minY]]
    271         trenchPts=numpy.array([minX+50., minY+50., -1000.])
     284        trenchPoly = [[minX+40., minY], [minX+40., minY+100.],
     285            [minX+60., minY+100.], [minX+60., minY]]
     286        trenchPts = numpy.array([minX+50., minY+50., -1000.])
    272287        #
    273         PtPolData=[[trenchPoly, trenchPts]]
    274         F=qs.quantity_from_Pt_Pol_Data_and_Raster(PtPolData, 'PointData_ElevTest.tif', domain)
    275 
    276         testPts_X=numpy.array([50., 3.])
    277         testPts_Y=numpy.array([1., 20.])
    278         fitted=F(testPts_X,testPts_Y)
     288        PtPolData = [[trenchPoly, trenchPts]]
     289        F = qs.quantity_from_Pt_Pol_Data_and_Raster(PtPolData,
     290            'PointData_ElevTest.tif', domain)
     291
     292        testPts_X = numpy.array([50., 3.])
     293        testPts_Y = numpy.array([1., 20.])
     294        fitted = F(testPts_X,testPts_Y)
     295
     296        import pdb
     297        pdb.set_trace()
    279298   
    280299        # The fitted value in the trench should be -1000.
    281         assert(fitted[0]==-1000.)
     300        assert(numpy.allclose(fitted[0],-1000.))
    282301
    283302        # Find the nearest domain point to the second test point
    284303        # This will have been used in constructing the elevation raster
    285         nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
    286         nearest_x=domain.centroid_coordinates[nearest,0]
     304        nearest = ((domain.centroid_coordinates[:,0]-3.)**2 +\
     305                   (domain.centroid_coordinates[:,1]-20.)**2).argmin()
     306        nearest_x = domain.centroid_coordinates[nearest,0]
    287307        assert(numpy.allclose(fitted[1],-nearest_x/150.))
    288308
Note: See TracChangeset for help on using the changeset viewer.