Changeset 9338 for trunk/anuga_core/source/anuga
- Timestamp:
- Sep 18, 2014, 4:47:46 PM (11 years ago)
- Location:
- trunk/anuga_core/source/anuga/utilities
- Files:
-
- 2 edited
-
quantity_setting_functions.py (modified) (9 diffs)
-
spatialInputUtil.py (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py
r9320 r9338 8 8 import anuga.utilities.spatialInputUtil as su 9 9 10 def make_nearestNeighbour_quantity_function(quantity_xyValueIn, domain, threshold_distance = 9.0e+100, background_value = 9.0e+100): 10 def make_nearestNeighbour_quantity_function( 11 quantity_xyValueIn, 12 domain, 13 threshold_distance = 9.0e+100, 14 background_value = 9.0e+100): 11 15 """ 12 16 Function which makes another function, which can be used in set_quantity … … 14 18 Idea: For every point x,y in the domain, we want to set a quantity based on 15 19 the 'nearest-neighbours' from quantity_xyValue (a 3 column array with 16 x,y,quantity-value), UNLESS the distance from x,y to the nearest-neighbour is > threshold_distance. 17 In the latter case, we want to set the quantity value to 'background_value' 20 x,y,quantity-value), 21 UNLESS the distance from x,y to the nearest-neighbour is > 22 threshold_distance. 23 In the latter case, we want to set the quantity value to 24 'background_value' 18 25 19 26 We need a function f(x,y) to do that. This routine makes the 20 function, with the desired quantity_xyValue points, threshold_distance, and21 background_value27 function, with the desired quantity_xyValue points, 28 threshold_distance, and background_value 22 29 23 30 INPUTS: … … 26 33 defining the points used to set the new quantity values 27 34 domain -- The ANUGA domain 28 threshold_distance -- Points greater than this distance from their nearest quantity_xyValue point are set to background_value 35 threshold_distance -- Points greater than this distance from their 36 nearest quantity_xyValue point are set to background_value 29 37 background_value -- see 'threshold_distance' 30 38 … … 42 50 quantity_xyValue=copy.copy(quantity_xyValueIn.reshape((1,3))) 43 51 # Make a function which gives us the ROW-INDEX of the nearest xy point in quantity_xyValue 44 quantity_xy_interpolator=scipy.interpolate.NearestNDInterpolator(quantity_xyValue[:,0:2], scipy.arange(len(quantity_xyValue[:,2]))) 52 quantity_xy_interpolator=scipy.interpolate.NearestNDInterpolator( 53 quantity_xyValue[:,0:2], 54 scipy.arange(len(quantity_xyValue[:,2]))) 45 55 46 56 # … … 67 77 q_index=quantity_xy_interpolator(z) 68 78 # Next find indices with distance < threshold_distance 69 dist_lt_thresh=( (z[:,0]-quantity_xyValue[q_index,0])**2 + (z[:,1]-quantity_xyValue[q_index,1])**2 < threshold_distance**2) 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) 70 82 dist_lt_thresh=dist_lt_thresh.nonzero()[0] 71 83 quantity_output[dist_lt_thresh] = quantity_xyValue[q_index[dist_lt_thresh],2] … … 78 90 def composite_quantity_setting_function(poly_fun_pairs, domain): 79 91 """ 80 Make a 'composite function' to set quantities -- applies different functions inside different polygon regions. 92 Make a 'composite function' to set quantities -- applies different 93 functions inside different polygon regions. 81 94 82 95 poly_fun_pairs = [ [p0, f0], [p1, f1], ...] … … 84 97 Where: 85 98 86 fi is a function, or a constant, or the name of a 87 gdal-compatible rasterFile, or a numpy array with 3 columns; 88 89 pi is a polygon, or a polygon filename (shapefile or a csv format that anuga.read_polygon will read), 90 or None, or 'All' (or it can be 'Extent' in the case that fi is a 91 rasterFile name) 99 fi is a function, 100 or a constant, 101 or a '.txt' or '.csv' file with comma separated xyz data 102 and a single header row, 103 or the name of a gdal-compatible rasterFile 104 (not ending in .txt or .csv), 105 or a numpy array with 3 columns; 106 107 pi is a polygon, 108 or a polygon filename (shapefile or a csv format that 109 anuga.read_polygon will read), 110 or None ( equivalent to a polygon with zero area), 111 or 'All' (equivalent to a polygon covering everything) 112 or 'Extent' in the case that fi is a rasterFile name 113 (equivalent to a polygon with the same extent as the raster) 92 114 93 IMPORTANT: When polygons overlap, the first elements of the list are given priority. 115 IMPORTANT: When polygons overlap, the first elements of the list are 116 given priority. 94 117 The approach is: 95 First f0 is applied to all points in p0, and we record that these points have been 'set' 96 Next f1 is applied to all points in p1 which have not been 'set', and then we record those points as being 'set' 97 Next f2 is applied to all points in p2 which have not been 'set', and then we record those points as being 'set' 118 First f0 is applied to all points in p0, and we record 119 that these points have been 'set' 120 Next f1 is applied to all points in p1 which have not 121 been 'set', and then we record those points as being 'set' 122 Next f2 is applied to all points in p2 which have not 123 been 'set', and then we record those points as being 'set' 98 124 ... etc 99 125 100 126 INPUT: 101 poly_fun_pairs = [ [p0, f0], [p1, f1], ...] 102 103 where fi(x,y) is a function returning quantity values at points, or any of the special cases below 104 SPECIAL fi CASES: 105 fi = a constant in which case points in the polygon are set to that value, 106 fi = a string rasterFile name which can be passed to quantityRasterFun to make a function, 107 fi = a numpy array with 3 columns (x,y,Value) in which case nearest-neighbour interpolation is used on the points 108 109 110 pi are polygons where we want to use fi inside 111 SPECIAL pi CASES: 112 If pi is a filename ending in .shp or a csv format that anuga.read_polygon can read, we assume it contains a polygon we have to read 113 If any pi = 'All', then we assume that ALL unset points are set 114 using the function. This CAN ONLY happen in the last [fi,pi] pair where pi is 115 not None (since fi will be applied to all remaining points -- so anything else is probably an input mistake) 116 If any pi = None, then that [fi,pi] pair is skipped 117 If pi = 'Extent' and fi is the name of a raster file, then the 118 extent of the raster file is used to define the polygon 119 120 domain = ANUGA domain object 127 poly_fun_pairs = [ [p0, f0], [p1, f1], ...] 128 129 where fi(x,y) is a function returning quantity values at points, 130 or any of the special cases below 131 SPECIAL fi CASES: 132 fi = a constant in which case points in the polygon are 133 set to that value, 134 fi = a .txt or .csv file name containing x, y, z data, 135 with comma separators and a single header row 136 fi = a string rasterFile name (not ending in .txt or .csv) 137 which can be passed to quantityRasterFun to make a function, 138 fi = a numpy array with 3 columns (x,y,Value) in which case 139 nearest-neighbour interpolation is used on the points 140 141 142 pi are polygons where we want to use fi inside 143 SPECIAL pi CASES: 144 If pi is a filename ending in .shp or a csv format that 145 anuga.read_polygon can read, we assume it contains a polygon we have to read 146 If any pi = 'All', then we assume that ALL unset points are set 147 using the function. This CAN ONLY happen in the last [fi,pi] pair where pi is 148 not None (since fi will be applied to all remaining points 149 -- so anything else is probably an input mistake) 150 If any pi = None, then that [fi,pi] pair is skipped 151 If pi = 'Extent' and fi is the name of a raster file, then the 152 extent of the raster file is used to define the polygon 153 154 domain = ANUGA domain object 121 155 122 156 … … 149 183 if (poly_fun_pairs[i][0]=='All'): 150 184 # This is only ok if all the othe poly_fun_pairs are None 151 remaining_poly_fun_pairs_are_None=[ poly_fun_pairs[j][0] is None for j in range(i+1,lfp)] 185 remaining_poly_fun_pairs_are_None = \ 186 [ poly_fun_pairs[j][0] is None for j in range(i+1,lfp)] 152 187 if(not all(remaining_poly_fun_pairs_are_None)): 153 188 raise Exception, 'Can only have the last polygon = All' … … 191 226 elif ( type(fi) is str and os.path.exists(fi)): 192 227 # fi is a file which is assumed to be 193 # a gdal-compatible raster 194 newfi = quantityRasterFun(domain, fi) 195 quantityVal[fInds] = newfi(x[fInds], y[fInds]) 228 # a gdal-compatible raster OR an x,y,z elevation file 229 if os.path.splitext(fi)[1] in ['.txt', '.csv']: 230 # Treating input file ' + fi + ' as xyz array with 1 header row 231 fi_array = numpy.genfromtxt(fi, delimiter=",", skip_header=1) 232 if fi_array.shape[1] is not 3: 233 print 'Treated input file ' + fi + ' as xyz array with 1 header row' 234 raise Exception, 'Array should have 3 columns -- x,y,value' 235 newfi = make_nearestNeighbour_quantity_function(fi_array, domain) 236 quantityVal[fInds] = newfi(x[fInds], y[fInds]) 237 else: 238 # Treating input file as a raster 239 newfi = quantityRasterFun(domain, fi) 240 quantityVal[fInds] = newfi(x[fInds], y[fInds]) 196 241 elif(type(fi) is numpy.ndarray): 197 242 if fi.shape[1] is not 3: -
trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py
r9335 r9338 251 251 except: 252 252 msg = 'Could not read points from ' + filename +\ 253 '. Make sure it has a single header row, with comma separator, and the first 2 columns are x,y' 253 '. Make sure it has a single header row, ' +\ 254 'with comma separator, and the first 2 columns are x,y' 254 255 raise Exception, msg 255 256 return points … … 408 409 as joined and so defining line segments 409 410 410 OUTPUT: The squared distance, and the index i of the segment [x_i,y_i],[x_i+1,y_i+1] in segments which is closest to pt 411 OUTPUT: The squared distance, and the index i of the segment 412 [x_i,y_i],[x_i+1,y_i+1] in segments which is closest to pt 411 413 """ 412 414 ll=len(segments) … … 425 427 def shift_point_on_line(pt, lineIn, nearest_segment_index): 426 428 """ 427 Support pt is a point, which is near to the 'nearest_segment_index' segment of the line 'lineIn' 429 Support pt is a point, which is near to the 'nearest_segment_index' 430 segment of the line 'lineIn' 428 431 429 432 This routine moves the nearest end point of that segment on line to pt. … … 431 434 INPUTS: pt -- [x,y] point 432 435 lineIn -- [ [x0, y0], [x1, y1], ..., [xN,yN]] 433 nearest_segment_index = index where the distance of pt to the line from [x_i,y_i] to [x_i+1,y_i+1] is minimum 436 nearest_segment_index = index where the distance of pt to 437 the line from [x_i,y_i] to [x_i+1,y_i+1] is minimum 434 438 435 439 OUTPUT: The new line … … 473 477 p0 = L1_pts[tmp[1]] 474 478 p1 = L1_pts[tmp[1]+1] 475 endPt_Dist_Sq=min( ( (p0[0]-iP[0])**2 + (p0[1]-iP[1])**2), ( (p1[0]-iP[0])**2 + (p1[1]-iP[1])**2)) 479 endPt_Dist_Sq=min( ( (p0[0]-iP[0])**2 + (p0[1]-iP[1])**2), 480 ( (p1[0]-iP[0])**2 + (p1[1]-iP[1])**2)) 476 481 # 477 482 if(endPt_Dist_Sq>point_movement_threshold**2): … … 560 565 tol2 = [see check_polygon_is_small] Probably doesn't need to change 561 566 562 nameFlag = This will be printed if intersection occurs. Convenient way to display intersecting filenames 567 nameFlag = This will be printed if intersection occurs. 568 Convenient way to display intersecting filenames 563 569 564 570 OUTPUTS: L1,L2 with intersection points added in the right places … … 567 573 if(L1.Intersects(L2)): 568 574 # Get points on the lines 569 L1_pts=Wkb2ListPts(L1) #L1.GetPoints()570 L2_pts=Wkb2ListPts(L2) #L2.GetPoints()575 L1_pts=Wkb2ListPts(L1) 576 L2_pts=Wkb2ListPts(L2) 571 577 572 578 # Buffer lines by a small amount … … 578 584 if(L1_L2_intersect.GetGeometryCount()==1): 579 585 if(not check_polygon_is_small(L1_L2_intersect, buf, tol2)): 580 #print L1_L2_intersect.GetEnvelope() 581 raise Exception, 'line intersection is not allowed. Envelope %s '% str(L1_L2_intersect.GetEnvelope()) 586 msg = 'line intersection is not allowed. ' + \ 587 'Envelope %s '% str(L1_L2_intersect.GetEnvelope()) 588 raise Exception(msg) 582 589 # Seems to need special treatment with only 1 intersection point 583 590 intersectionPts=[L1_L2_intersect.Centroid().GetPoint()] … … 597 604 # Insert the points into the line segments 598 605 for i in range(len(intersectionPts)): 599 L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, point_movement_threshold, verbose=verbose) 600 L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, point_movement_threshold, verbose=verbose) 606 L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, 607 point_movement_threshold, verbose=verbose) 608 L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, 609 point_movement_threshold, verbose=verbose) 601 610 602 611 # Convert to the input format … … 716 725 get a range of raster values inside a polygon 717 726 718 Approach: A 'trial-grid' of points is created which is 'almost' covering the range of the polygon (xmin-xmax,ymin-ymax). 727 Approach: A 'trial-grid' of points is created which is 'almost' 728 covering the range of the polygon (xmin-xmax,ymin-ymax). 719 729 720 730 (Actually it is just inside this region, to avoid polygon-boundary issues, see below) … … 763 773 keepers = inside_polygon(Grid, polygon) 764 774 if(len(keepers)==0): 765 import pdb 766 pdb.set_trace() 767 ## Now find the xy values which are actually inside the polygon 768 #polpath=path.Path(polygonArr) 769 #keepers=[] 770 #for i in range(len(Grid[:,0])): 771 # if(polpath.contains_point(Grid[i,:])): 772 # keepers=keepers+[i] 773 ##import pdb 774 ##pdb.set_trace() 775 raise Exception('No points extracted from polygon') 775 776 xyInside=Grid[keepers,:] 776 777 … … 958 959 msg='Breakline/riverwall point '+str(blCat[n1][j][0:2]) +' on '+ n1+\ 959 960 ' is outside the bounding polygon.\n'+\ 960 'Check that it exceeds the bounding polygon by a distance < point_movement_threshold \n'+\ 961 'Check that it exceeds the bounding polygon'+\ 962 ' by a distance < point_movement_threshold \n'+\ 961 963 ' so it can be moved back onto the polygon' 962 964 print 'Polygon\n ' … … 1013 1015 Read them in 1014 1016 1015 INPUT: fileList -- a list of shapefile and/or anuga_polygon csv filenames [e.g. from glob.glob('GIS/Breaklines/*.shp')] 1017 INPUT: fileList -- a list of shapefile and/or anuga_polygon csv filenames 1018 [e.g. from glob.glob('GIS/Breaklines/*.shp')] 1016 1019 1017 1020 OUTPUT: dictionary with breaklines [filenames are keys] … … 1027 1030 def readListOfRiverWalls(rwfileList): 1028 1031 """ 1029 Take a list with the names of riverwall input files [should be comma-separated x,y,elevation files] 1030 1031 The input file can OPTIONALLY have a first line defining some hydraulic parameters. A valid example 1032 is 1032 Take a list with the names of riverwall input files 1033 [should be comma-separated x,y,elevation files] 1034 1035 The input file can OPTIONALLY have a first line defining some 1036 hydraulic parameters. A valid example is 1033 1037 1034 1038 Qfactor: 1.5, s1: 0.94 … … 1039 1043 Read their coordinates into a dict with their names, read for use by ANUGA 1040 1044 1041 INPUT: rwfileList -- a list of riverwall filenames [e.g. from glob.glob('GIS/RiverWalls/*.csv')] 1045 INPUT: rwfileList -- a list of riverwall filenames 1046 [e.g. from glob.glob('GIS/RiverWalls/*.csv')] 1042 1047 1043 1048 OUTPUT: … … 1074 1079 def combine_breakLines_and_riverWalls_for_mesh(breakLines, riverWalls): 1075 1080 """ 1076 Combine breaklines and riverwalls for input in mesh generation, ensuring both have 2 coordinates only 1081 Combine breaklines and riverwalls for input in mesh generation, 1082 ensuring both have 2 coordinates only 1077 1083 """ 1078 1084 mesh_breakLines=riverWalls.values()+breakLines.values() 1079 1085 for i in range(len(mesh_breakLines)): 1080 mesh_breakLines[i] = [mesh_breakLines[i][j][0:2] for j in range(len(mesh_breakLines[i]))] 1086 mesh_breakLines[i] =\ 1087 [mesh_breakLines[i][j][0:2] for j in range(len(mesh_breakLines[i]))] 1081 1088 return mesh_breakLines 1082 1089 … … 1089 1096 Can do this with the current function 1090 1097 1091 INPUTS: pattern == character string containing pattern which is inside exactly 2 keys in breaklines 1098 INPUTS: pattern == character string containing pattern which 1099 is inside exactly 2 keys in breaklines 1092 1100 1093 1101 breakLines = the breakLinesIn dictionary 1094 1102 1095 reverse2nd = True/False or None. Reverse the order of the 2nd set of edges before making the polygon. 1103 reverse2nd = True/False or None. Reverse the order of the 1104 2nd set of edges before making the polygon. 1096 1105 If None, then we compute the distance between the 1097 first point on breakline1 and the first/last points on breakline2, and reverse2nd if 1098 the 'distance from the first point' < 'distance from the last point' 1106 first point on breakline1 and the first/last 1107 points on breakline2, and reverse2nd if 1108 the 'distance from the first point' < 1109 'distance from the last point' 1099 1110 1100 1111 OUTPUT: Polygon made with the 2 breaklines … … 1129 1140 breakLines[bk[matchers[1]]].reverse() 1130 1141 polyOut=breakLines[bk[matchers[0]]] + breakLines[bk[matchers[1]]] 1131 #polyOut=polyOut+[polyOut[0]] 1132 1133 # If the breakLines values have > 2 entries (i.e. like for riverwalls),remove the third1142 1143 # If the breakLines values have > 2 entries (i.e. like for riverwalls), 1144 # remove the third 1134 1145 if(len(polyOut[0])>2): 1135 1146 polyOut=[polyOut[i][0:2] for i in range(len(polyOut))]
Note: See TracChangeset
for help on using the changeset viewer.
