Changeset 9513
- Timestamp:
- Jan 27, 2015, 5:03:44 PM (10 years ago)
- Location:
- trunk/anuga_core/source/anuga/utilities
- Files:
-
- 4 edited
-
quantity_setting_functions.py (modified) (20 diffs)
-
spatialInputUtil.py (modified) (5 diffs)
-
test/test_quantity_setting_functions.py (modified) (3 diffs)
-
test/test_spatialInputUtil.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py
r9511 r9513 32 32 INPUTS: 33 33 @param quantity_xyValueIn -- A 3 column array with 'x,y, Value' 34 defining the points used to set the new quantity values 34 defining the points used to set the new quantity values in 35 georeferenced coordinates 35 36 @param domain -- The ANUGA domain 36 37 @param k_nearest_neighbors If > 1, then an inverse-distance-weighted average … … 131 132 domain, 132 133 clip_range = None, 133 nan_treatment = 'exception'): 134 nan_treatment = 'exception', 135 nan_interpolation_region_polygon = None, 136 default_k_nearest_neighbours = 1, 137 default_raster_interpolation = 'pixel', 138 verbose=True): 134 139 """ Make a 'composite function' to set quantities -- applies different 135 140 functions inside different polygon regions. … … 176 181 fi = a .txt or .csv file name containing x, y, z data, 177 182 with comma separators and an optional header row 178 containing letters 183 containing letters (nearest neighbour interpolation is used) 179 184 fi = a string rasterFile name (not ending in .txt or .csv) 180 which can be passed to quantityRasterFun to make a function ,185 which can be passed to quantityRasterFun to make a function 181 186 fi = a numpy array with 3 columns (x,y,Value) in which case 182 187 nearest-neighbour interpolation is used on the points … … 212 217 poly,fun pairs (in sequence) to set the value. 213 218 219 @param nan_interpolation_region_polygon = None, or 'All', or a list 220 of csv or shp filenames containing polygons, or a list of 221 anuga polygon objects. 222 223 If it is not None, then all x,y points which evaluate to nan 224 on their **first preference** dataset are recorded, and as a 225 final step, the values at these x,y points 226 **which are inside the nan_interpolation_region_polygon** 227 are interpolated from the other x,y,F(x,y) values. 228 229 Nearest neighbour interpolation is used, with 230 k_nearest_neighbours taken from default_k_nearest_neighbours. 231 232 Note that if nan_treatment = 'exception', then nan's will cause 233 exceptions earlier on in this routine, so you will need 234 nan_treatment = 'fall_through' to use this option. 235 236 Example of why you might want this: 237 Say you have 2 elevation datasets (one defining the 238 topography above MSL, and the other defining the topography 239 below MSL). There might be small nan gaps between them, 240 which you would like to fill with interpolation. That 241 can be done with this option. 242 243 @param default_k_nearest_neighbours = integer >=1 . The value of 244 k_nearest_neighbours passed to 245 make_nearestNeighbour_quantity_function when a 'special_case' 246 value of fi is passed in (either a point array or a .txt or 247 .csv point file), or when nan_interpolation_region_polygon is 248 not None 249 250 @param default_raster_interpolation = 'pixel' or 'bilinear'. The value of 251 'interpolation' passed to quantityRasterFun if a raster filename 252 is passed as one of the fi. 253 254 @param verbose TRUE/FALSE Print more information 255 214 256 OUTPUT: A function F(x,y) which can be used e.g. to set the quantity 215 257 domain.set_quantity('elevation', F) 258 216 259 """ 217 260 import os … … 230 273 raise Exception('clip_range minima must be less than maxima') 231 274 275 232 276 def F(x,y): 233 """This is the function we return277 """This is the function returned by composite_quantity_setting_function 234 278 It can be passed to set_quantity 235 279 """ 236 isSet=numpy.zeros(len(x)) # 0/1 - record if each point has been set 237 quantityVal=x*0 # Function return value 238 lpf=len(poly_fun_pairs) 239 if(lpf<=0): 280 isSet = numpy.zeros(len(x)) # 0/1 - record if each point has been set 281 quantityVal = x*0 + numpy.nan # Function return value 282 283 # Record points which evaluated to nan on their first preference 284 # dataset. 285 was_ever_nan = (x*0).astype(int) 286 287 lpf = len(poly_fun_pairs) 288 if(lpf <= 0): 240 289 raise Exception('Must have at least 1 fun-poly-pair') 241 290 242 291 # Make an array of 'transformed' spatial coordinates, for checking 243 292 # polygon inclusion 244 xll =domain.geo_reference.xllcorner245 yll =domain.geo_reference.yllcorner246 xy_array_trans =numpy.vstack([x+xll,y+yll]).transpose()293 xll = domain.geo_reference.xllcorner 294 yll = domain.geo_reference.yllcorner 295 xy_array_trans = numpy.vstack([x+xll,y+yll]).transpose() 247 296 248 297 # Check that none of the pi polygons [except perhaps the last] is 'All' … … 255 304 raise Exception('Can only have the last polygon = All') 256 305 257 # M AIN LOOP306 # Main Loop 258 307 # Apply the fi inside the pi 259 308 for i in range(lpf): … … 271 320 if(pi == 'All'): 272 321 # Get all unset points 273 fInside =(1-isSet)274 fInds =(fInside==1).nonzero()[0]322 fInside = (1-isSet) 323 fInds = (fInside==1).nonzero()[0] 275 324 276 325 else: … … 289 338 290 339 # Then we get the extent from the raster itself 291 pi_path =su.getRasterExtent(fi,asPolygon=True)340 pi_path = su.getRasterExtent(fi,asPolygon=True) 292 341 293 342 elif(type(pi)==str and os.path.isfile(pi) ): 294 343 # pi is a file 295 pi_path =su.read_polygon(pi)344 pi_path = su.read_polygon(pi) 296 345 297 346 else: 298 347 # pi is the actual polygon data 299 pi_path =pi348 pi_path = pi 300 349 301 350 # Get the insides of unset points inside pi_path 302 notSet =(isSet==0.).nonzero()[0]351 notSet = (isSet==0.).nonzero()[0] 303 352 fInds = inside_polygon(xy_array_trans[notSet,:], pi_path) 304 353 fInds = notSet[fInds] … … 320 369 elif isinstance(fi, (int, long, float)): 321 370 # fi is a numerical constant 322 quantityVal[fInds] =fi*1.0371 quantityVal[fInds] = fi*1.0 323 372 324 373 elif ( type(fi) is str and os.path.exists(fi)): … … 330 379 if fi_array.shape[1] is not 3: 331 380 print 'Treated input file ' + fi +\ 332 ' as xyz array with '+ \ 333 str(int(hasLetters)) + ' header row' 381 ' as xyz array with an optional header' 334 382 msg = 'Array should have 3 columns -- x,y,value' 335 383 raise Exception(msg) 336 384 337 385 newfi = make_nearestNeighbour_quantity_function( 338 fi_array, domain) 386 fi_array, domain, 387 k_nearest_neighbours = default_k_nearest_neighbours) 339 388 quantityVal[fInds] = newfi(x[fInds], y[fInds]) 340 389 341 390 else: 342 391 # Treating input file as a raster 343 newfi = quantityRasterFun(domain, fi) 392 newfi = quantityRasterFun(domain, fi, 393 interpolation = default_raster_interpolation) 344 394 quantityVal[fInds] = newfi(x[fInds], y[fInds]) 345 395 … … 348 398 msg = 'Array should have 3 columns -- x,y,value' 349 399 raise Exception(msg) 350 newfi = make_nearestNeighbour_quantity_function(fi, domain) 400 newfi = make_nearestNeighbour_quantity_function(fi, domain, 401 k_nearest_neighbours = default_k_nearest_neighbours) 351 402 quantityVal[fInds] = newfi(x[fInds], y[fInds]) 352 403 … … 354 405 print 'Error with function from' 355 406 print fi 356 msg='Cannot make function from type ' +str(type(fi))407 msg='Cannot make function from type ' + str(type(fi)) 357 408 raise Exception, msg 358 409 … … 361 412 ################################################################### 362 413 363 nan_flag = quantityVal[fInds] != quantityVal[fInds]414 nan_flag = (quantityVal[fInds] != quantityVal[fInds]) 364 415 nan_inds = nan_flag.nonzero()[0] 416 was_ever_nan[fInds[nan_inds]] = 1 417 365 418 if len(nan_inds)>0: 366 419 if nan_treatment == 'exception': … … 379 432 'in composite_quantity_setting_function. ' + \ 380 433 'They will be passed to later poly_fun_pairs' 381 print msg434 if verbose: print msg 382 435 not_nan_inds = (1-nan_flag).nonzero()[0] 383 436 … … 388 441 msg = '( Actually all the values were nan - ' + \ 389 442 'Are you sure they should be? Possible error?)' 390 print msg443 if verbose: print msg 391 444 continue 445 392 446 else: 393 447 msg = 'Found nan values in ' + \ … … 397 451 398 452 # Record that the points have been set 399 isSet[fInds] =1453 isSet[fInds] = 1 400 454 401 455 # Enforce clip_range … … 410 464 # End of loop 411 465 466 # Find points which were nan on their first preference dataset + are 467 # inside nan_interpolation_region_polygon. Then reinterpolate their 468 # values from the other x,y, quantityVal points. 469 if (nan_interpolation_region_polygon is not None) &\ 470 (was_ever_nan.sum() > 0): 471 if nan_interpolation_region_polygon == 'All': 472 points_to_reinterpolate = was_ever_nan.nonzero()[0] 473 else: 474 # nan_interpolation_region_polygon contains information on 1 or 475 # more polygons 476 # Inside those polygons, we need to re-interpolate points which 477 # first evaluted to na 478 possible_points_to_reint = was_ever_nan.nonzero()[0] 479 points_to_reinterpolate = numpy.array([]).astype(int) 480 481 for i in range(len(nan_interpolation_region_polygon)): 482 nan_pi = nan_interpolation_region_polygon[i] 483 484 # Ensure nan_pi = list of x,y points making a polygon 485 if(type(nan_pi) == str): 486 nan_pi = su.read_polygon(nan_pi) 487 488 points_in_nan_pi = inside_polygon( 489 xy_array_trans[possible_points_to_reint,:], 490 nan_pi) 491 492 if len(points_in_nan_pi)>0: 493 points_to_reinterpolate = numpy.hstack( 494 [points_to_reinterpolate, 495 points_in_nan_pi]) 496 497 if verbose: 498 print 'Re-interpolating ', len(points_to_reinterpolate),\ 499 ' points which were nan under their',\ 500 ' first-preference and are inside the',\ 501 ' nan_interpolation_region_polygon' 502 503 # Find the interpolation points = points not needing reinterpolation 504 ip = x*0 + 1 505 ip[points_to_reinterpolate] = 0 506 number_of_ip = ip.sum() 507 ip = ip.nonzero()[0] 508 509 if(number_of_ip < default_k_nearest_neighbours): 510 raise Exception('Too few non-nan points to interpolate from') 511 512 # Make function for re-interpolation. Note this requires 513 # x,y,z in georeferenced coordinates, whereas x,y are ANUGA 514 # coordinates 515 reinterp_F = make_nearestNeighbour_quantity_function( 516 numpy.vstack([xy_array_trans[ip,0], xy_array_trans[ip,1], 517 quantityVal[ip]]).transpose(), 518 domain, 519 k_nearest_neighbours = default_k_nearest_neighbours) 520 521 # re-interpolate 522 quantityVal[points_to_reinterpolate] = reinterp_F( 523 x[points_to_reinterpolate], y[points_to_reinterpolate]) 524 525 isSet[points_to_reinterpolate] = 1 526 527 # Check there are no remaining nan values 412 528 if( min(isSet) != 1): 413 print 'Some points were not inside any polygon, ',\ 414 'or evalute to nan over all datasets' 529 print 'Some points remain as nan, which is not allowed' 415 530 unset_inds = (isSet!=1).nonzero()[0] 416 531 lui = min(5, len(unset_inds)) … … 428 543 ############################################################################## 429 544 430 def quantityRasterFun(domain, rasterFile ):545 def quantityRasterFun(domain, rasterFile, interpolation='pixel'): 431 546 """ 432 547 Make a function whick takes x,y in ANUGA coordinates, and returns the values 433 548 on a raster rasterFile 434 549 435 This can be used to set a quantity, and takes care of the manual conversion from 436 ANUGA coordinates to spatial coordinates. 437 438 INPUTS: domain = ANUGA domain 439 rasterFile = Filename of the raster to extract point values from 550 This can be used to set a quantity, and takes care of the manual conversion 551 from ANUGA coordinates to spatial coordinates. 552 553 INPUTS: @param domain = ANUGA domain 554 @param rasterFile = Filename of the raster to extract point values 555 from 556 @param interpolation = 'pixel' (in which case the point value is 557 set from the pixel it is on) or 'bilinear' in which case 558 the point value is set from bilinear interpolation of 559 pixels. 440 560 441 561 OUTPUT: Function which takes x,y in ANUGA coordinates, and outputs their … … 448 568 yll=domain.geo_reference.yllcorner 449 569 inDat=scipy.vstack([x+xll,y+yll]).transpose() 450 return rasterValuesAtPoints(xy=inDat,rasterFile=rasterFile) 570 return rasterValuesAtPoints(xy=inDat,rasterFile=rasterFile, 571 interpolation=interpolation) 451 572 452 573 return QFun … … 461 582 interpolator is used. 462 583 584 This has been superceeded by composite_quantity_setting_function 585 463 586 INPUT: 464 Pt_Pol_Data = a list with [ [ Polygon_0, Pt_XYZ_0], 465 [ Polygon_1, Pt_XYZ_1], 466 ... ] 587 @param Pt_Pol_Data = a list with [ [ Polygon_0, Pt_XYZ_0], 588 [ Polygon_1, Pt_XYZ_1], 589 ... 590 ] 467 591 Here Polygon_i is a polygon in ANUGA format, 468 592 and Pt_XYZ_i is a 3 column array of x,y,Value points 469 quantity_raster = A GDAL-compatible quantity raster470 domain = ANUGA domain593 @param quantity_raster = A GDAL-compatible quantity raster 594 @param domain = ANUGA domain 471 595 """ 472 596 -
trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py
r9373 r9513 696 696 rasterFile, 697 697 band=1, 698 nodata_rel_tol = 1.0e-08): 698 nodata_rel_tol = 1.0e-08, 699 interpolation = 'pixel'): 699 700 """ 700 701 Get raster values at point locations. … … 712 713 This allows for truncation errors in nodata values which seem 713 714 to be introduced by some file-type conversions 715 716 @param interpolation 'pixel' or 'bilinear' determines how the 717 raster cell values are used to set the point value 714 718 715 719 OUTPUT: … … 730 734 731 735 # Get coordinates in pixel values 732 px = ( (xy[:,0] - xOrigin) / pixelWidth).astype(int) #x733 py = ( (xy[:,1] - yOrigin) / pixelHeight).astype(int) #y736 px = (xy[:,0] - xOrigin) / pixelWidth 737 py = (xy[:,1] - yOrigin) / pixelHeight 734 738 735 739 # Hold elevation 736 elev =px*0.740 elev = px*0. 737 741 738 742 # Get the right character for struct.unpack … … 753 757 754 758 # Upper bounds for pixel values, so we can fail gracefully 755 xMax =raster.RasterXSize756 yMax =raster.RasterYSize759 xMax = raster.RasterXSize 760 yMax = raster.RasterYSize 757 761 if(px.max()<xMax and px.min()>=0 and py.max()<yMax and py.min()>=0): 758 762 pass … … 763 767 # Get values -- seems we have to loop, but it is efficient enough 764 768 for i in range(len(px)): 765 xC=int(px[i]) 766 yC=int(py[i]) 767 structval=rasterBand.ReadRaster(xC,yC,1,1,buf_type=rasterBand.DataType) 768 elev[i] = struct.unpack(CtypeName, structval)[0] 769 770 if(interpolation=='pixel'): 771 # Pixel coordinates 772 xC=int(numpy.floor(px[i])) 773 yC=int(numpy.floor(py[i])) 774 775 structval = rasterBand.ReadRaster(xC,yC,1,1, 776 buf_type=rasterBand.DataType) 777 elev[i] = struct.unpack(CtypeName, structval)[0] 778 779 elif(interpolation=='bilinear'): 780 # Pixel coordinates 781 xl = int(numpy.floor(px[i])) 782 yl = int(numpy.floor(py[i])) 783 784 # Find neighbours required for bilinear interpolation 785 # l = lower, u = upper 786 if(px[i] - xl > 0.5): 787 xu = min(xl + 1, xMax - 1) 788 else: 789 # Swap xl for xu 790 xu = xl 791 xl = max(xu - 1, 0) 792 793 if(py[i] - yl > 0.5): 794 yu = min(yl + 1, yMax - 1) 795 else: 796 yu = yl 797 yl = max(yu - 1, 0) 798 799 # Map x,y to unit square 800 x = px[i] - (xl + 0.5) 801 y = py[i] - (yl + 0.5) 802 803 if not ( (x>=0.) & (x<=1.)): 804 print x, xl, xu, px[i] 805 raise Exception('x out of bounds') 806 807 if not ( (y>=0.) & (y<=1.)): 808 print y, yl, yu, py[i] 809 raise Exception('y out of bounds') 810 811 # Lower-left 812 structval = rasterBand.ReadRaster(xl,yl,1,1, 813 buf_type=rasterBand.DataType) 814 r00 = struct.unpack(CtypeName, structval)[0] 815 # Upper left 816 structval = rasterBand.ReadRaster(xl,yu,1,1, 817 buf_type=rasterBand.DataType) 818 r01 = struct.unpack(CtypeName, structval)[0] 819 # Lower-right 820 structval = rasterBand.ReadRaster(xu,yl,1,1, 821 buf_type=rasterBand.DataType) 822 r10 = struct.unpack(CtypeName, structval)[0] 823 # Upper right 824 structval = rasterBand.ReadRaster(xu,yu,1,1, 825 buf_type=rasterBand.DataType) 826 r11 = struct.unpack(CtypeName, structval)[0] 827 828 # Bilinear interpolation 829 elev[i] = r00*(1-x)*(1-y) + r01*(1-x)*y +\ 830 r10*x*(1-y) + r11*x*y 831 else: 832 raise Exception('Unknown value of "interpolation"') 769 833 770 834 # Deal with nodata -
trunk/anuga_core/source/anuga/utilities/test/test_quantity_setting_functions.py
r9512 r9513 109 109 110 110 # Test that F evaluated in 'ANUGA coordinates' [lower-left = 0,0] is correct 111 xGet=numpy.array([0., 10., 10., 0., 100.]) + 0.1 112 yGet=numpy.array([0., 0., 20.,90., 100. ]) + 0.1 113 expected=numpy.floor(xGet+yGet) 114 output=F(xGet,yGet) 115 assert(numpy.allclose(output,expected)) 116 117 118 # Test with 3 nearest neighbours at points which are exactly on data points 119 F = qs.make_nearestNeighbour_quantity_function(inPts, domain, 120 threshold_distance = 9.0e+100, background_value = 9.0e+100, 121 k_nearest_neighbours = 3) 122 111 123 xGet=numpy.array([0., 10., 10., 0., 100.]) 112 124 yGet=numpy.array([0., 0., 20.,90., 100. ]) … … 114 126 output=F(xGet,yGet) 115 127 assert(numpy.allclose(output,expected)) 128 129 # Test with 4 nearest neighbours, at points which are not exactly on data points 130 F = qs.make_nearestNeighbour_quantity_function(inPts, domain, 131 threshold_distance = 9.0e+100, background_value = 9.0e+100, 132 k_nearest_neighbours = 4) 133 # Test at non-trivial location 134 xGet=numpy.array([0., 10., 10., 0., 99.]) + 0.5 135 yGet=numpy.array([0., 0., 20.,90., 99. ]) + 0.5 136 expected=xGet+yGet 137 output=F(xGet,yGet) 138 assert(numpy.allclose(output,expected)) 139 116 140 117 141 ################################################################################# … … 259 283 self.assertRaises(Exception, lambda: should_fail_1()) 260 284 285 ########################################################################### 286 # This example features a function with some nan return values, and uses 287 # the nan_interpolation_region_polygon to try to fix it 288 def f0(x,y): 289 output = x/10. 290 output[0] = numpy.nan 291 return output 292 293 F = qs.composite_quantity_setting_function( 294 [[trenchPoly, f0], ['Extent', 'PointData_ElevTest.tif']], 295 domain, 296 nan_treatment = 'fall_through', 297 nan_interpolation_region_polygon = [trenchPoly], 298 default_k_nearest_neighbours = 3, 299 default_raster_interpolation = 'bilinear', 300 verbose=False) 301 302 # Points where we test the function. We deliberately use many points with x=50, 303 # which happens to ensure that the nan value is replaced with the same 304 # value it would have had anyway 305 testPts_X = numpy.array([50.,50.00, 50., 50., 97., 51., 3.]) 306 testPts_Y = numpy.array([1., 2., 3. , 4 , 20., 30., 60.]) 307 fitted = F(testPts_X,testPts_Y) 308 309 # We should have no nan values 310 assert(sum(fitted!=fitted) == 0) 311 312 # Now the fitted value in the trench should be determined by f0 because 313 # the re-interpolation of nan values was designed to ensure it 314 assert(numpy.allclose(fitted[0],50./10.)) 315 261 316 return 262 317 -
trunk/anuga_core/source/anuga/utilities/test/test_spatialInputUtil.py
r9490 r9513 340 340 z_predicted=numpy.round(xA)+numpy.round(yA)-tifRange[0]-tifRange[2]-1.0 341 341 InDat=numpy.vstack([xA,yA]).transpose() 342 342 343 z_fitted=su.rasterValuesAtPoints(InDat, rasterFile='PointData_TestData.tif') 343 344 try: … … 345 346 except: 346 347 raise Exception, 'Error could be in rasterValuesAtPoints or in Make_Geotif' 348 349 350 # Try with bilinear interpolation 351 z_fitted=su.rasterValuesAtPoints(InDat, rasterFile='PointData_TestData.tif', 352 interpolation='bilinear') 353 z_predicted = xA + yA - tifRange[0] - tifRange[2] - 1.0 354 try: 355 assert(numpy.allclose(z_fitted,z_predicted)) 356 except: 357 raise Exception, 'Error could be in rasterValuesAtPoints or in Make_Geotif' 358 359 347 360 return 348 361
Note: See TracChangeset
for help on using the changeset viewer.
