Changeset 1911
- Timestamp:
- Oct 14, 2005, 6:44:45 AM (18 years ago)
- Location:
- inundation
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pyvolution/combine_pts.py
r1751 r1911 8 8 Geoscience Australia, 2005. 9 9 """ 10 from util import outside_polygon, inside_polygon10 from utilities.polygon import outside_polygon, inside_polygon 11 11 from Numeric import take, concatenate 12 12 import time -
inundation/pyvolution/data_manager.py
r1875 r1911 1261 1261 1262 1262 from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 1263 1263 from utilities.polygon import inside_polygon 1264 1264 1265 msg = 'Format must be either asc or ers' 1265 1266 assert format.lower() in ['asc', 'ers'], msg … … 1399 1400 #Interpolate 1400 1401 from least_squares import Interpolation 1401 from util import inside_polygon1402 1402 1403 1403 1404 #FIXME: This should be done with precrop = True (?), otherwise it'll … … 2880 2881 from Numeric import array, Float, concatenate, NewAxis, zeros,\ 2881 2882 sometrue 2882 2883 from utilities.polygon import inside_polygon 2883 2884 2884 2885 #FIXME: Should be variable … … 2984 2985 #Interpolate 2985 2986 from least_squares import Interpolation 2986 from util import inside_polygon2987 2987 2988 2988 2989 #FIXME: This should be done with precrop = True, otherwise it'll … … 3097 3098 from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 3098 3099 import ermapper_grids 3100 from utilities.polygon import inside_polygon 3099 3101 3100 3102 header = {} … … 3202 3204 #Interpolate 3203 3205 from least_squares import Interpolation 3204 from util import inside_polygon3206 3205 3207 3206 3208 #FIXME: This should be done with precrop = True (?), otherwise it'll -
inundation/pyvolution/least_squares.py
r1907 r1911 434 434 435 435 from quad import build_quadtree 436 from util import ensure_numeric 436 from utilities.numerical_tools import ensure_numeric 437 from utilities.polygon import inside_polygon 438 437 439 438 440 if data_origin is None: … … 489 491 if precrop is True: 490 492 from Numeric import take 491 from util import inside_polygon492 493 493 494 if verbose: print 'Getting boundary polygon' -
inundation/pyvolution/test_mesh.py
r1670 r1911 10 10 from config import epsilon 11 11 from Numeric import allclose, array 12 13 from utilities.polygon import inside_polygon 12 14 13 15 def distance(x, y): … … 705 707 from mesh import Mesh 706 708 from Numeric import zeros, Float 707 from util import inside_polygon708 709 709 710 #Create basic mesh … … 726 727 from mesh import Mesh 727 728 from Numeric import zeros, Float 728 from util import inside_polygon729 729 730 730 731 #Points … … 766 767 from mesh import Mesh 767 768 from Numeric import zeros, Float 768 from util import inside_polygon 769 769 770 770 771 #Points … … 807 808 from mesh import Mesh 808 809 from Numeric import zeros, Float 809 from util import inside_polygon810 810 from mesh_factory import rectangular 811 811 -
inundation/pyvolution/test_util.py
r1903 r1911 172 172 assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0) 173 173 174 def test_ensure_numeric(self):175 from util import ensure_numeric176 from Numeric import ArrayType, Float, array177 178 A = [1,2,3,4]179 B = ensure_numeric(A)180 assert type(B) == ArrayType181 assert B.typecode() == 'l'182 assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4183 184 185 A = [1,2,3.14,4]186 B = ensure_numeric(A)187 assert type(B) == ArrayType188 assert B.typecode() == 'd'189 assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4190 191 192 A = [1,2,3,4]193 B = ensure_numeric(A, Float)194 assert type(B) == ArrayType195 assert B.typecode() == 'd'196 assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0197 198 199 A = [1,2,3,4]200 B = ensure_numeric(A, Float)201 assert type(B) == ArrayType202 assert B.typecode() == 'd'203 assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0204 205 206 A = array([1,2,3,4])207 B = ensure_numeric(A)208 assert type(B) == ArrayType209 assert B.typecode() == 'l'210 assert A == B211 assert A is B #Same object212 213 214 A = array([1,2,3,4])215 B = ensure_numeric(A, Float)216 assert type(B) == ArrayType217 assert B.typecode() == 'd'218 assert A == B219 assert A is not B #Not the same object220 221 174 222 175 … … 1020 973 1021 974 1022 #Polygon stuff1023 def test_polygon_function_constants(self):1024 p1 = [[0,0], [10,0], [10,10], [0,10]]1025 p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]1026 1027 f = Polygon_function( [(p1, 1.0)] )1028 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p11029 assert allclose(z, [1,1,0,0])1030 1031 1032 f = Polygon_function( [(p2, 2.0)] )1033 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p21034 assert allclose(z, [2,0,0,2])1035 1036 1037 #Combined1038 f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )1039 z = f([5, 5, 27, 35], [5, 9, 8, -5])1040 assert allclose(z, [2,1,0,2])1041 1042 1043 def test_polygon_function_callable(self):1044 """Check that values passed into Polygon_function can be callable1045 themselves.1046 """1047 p1 = [[0,0], [10,0], [10,10], [0,10]]1048 p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]1049 1050 f = Polygon_function( [(p1, test_function)] )1051 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p11052 assert allclose(z, [10,14,0,0])1053 1054 #Combined1055 f = Polygon_function( [(p1, test_function), (p2, 2.0)] )1056 z = f([5, 5, 27, 35], [5, 9, 8, -5])1057 assert allclose(z, [2,14,0,2])1058 1059 1060 #Combined w default1061 f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)1062 z = f([5, 5, 27, 35], [5, 9, 8, -5])1063 assert allclose(z, [2,14,3.14,2])1064 1065 1066 #Combined w default func1067 f = Polygon_function( [(p1, test_function), (p2, 2.0)],1068 default = test_function)1069 z = f([5, 5, 27, 35], [5, 9, 8, -5])1070 assert allclose(z, [2,14,35,2])1071 1072 1073 def test_point_on_line(self):1074 1075 #Endpoints first1076 assert point_on_line( 0, 0, 0,0, 1,0 )1077 assert point_on_line( 1, 0, 0,0, 1,0 )1078 1079 #Then points on line1080 assert point_on_line( 0.5, 0, 0,0, 1,0 )1081 assert point_on_line( 0, 0.5, 0,1, 0,0 )1082 assert point_on_line( 1, 0.5, 1,1, 1,0 )1083 assert point_on_line( 0.5, 0.5, 0,0, 1,1 )1084 1085 #Then points not on line1086 assert not point_on_line( 0.5, 0, 0,1, 1,1 )1087 assert not point_on_line( 0, 0.5, 0,0, 1,1 )1088 1089 #From real example that failed1090 assert not point_on_line( 40,50, 40,20, 40,40 )1091 1092 1093 #From real example that failed1094 assert not point_on_line( 40,19, 40,20, 40,40 )1095 1096 1097 1098 1099 def test_inside_polygon_main(self):1100 1101 1102 #Simplest case: Polygon is the unit square1103 polygon = [[0,0], [1,0], [1,1], [0,1]]1104 1105 assert inside_polygon( (0.5, 0.5), polygon )1106 assert not inside_polygon( (0.5, 1.5), polygon )1107 assert not inside_polygon( (0.5, -0.5), polygon )1108 assert not inside_polygon( (-0.5, 0.5), polygon )1109 assert not inside_polygon( (1.5, 0.5), polygon )1110 1111 #Try point on borders1112 assert inside_polygon( (1., 0.5), polygon, closed=True)1113 assert inside_polygon( (0.5, 1), polygon, closed=True)1114 assert inside_polygon( (0., 0.5), polygon, closed=True)1115 assert inside_polygon( (0.5, 0.), polygon, closed=True)1116 1117 assert not inside_polygon( (0.5, 1), polygon, closed=False)1118 assert not inside_polygon( (0., 0.5), polygon, closed=False)1119 assert not inside_polygon( (0.5, 0.), polygon, closed=False)1120 assert not inside_polygon( (1., 0.5), polygon, closed=False)1121 1122 1123 1124 #From real example (that failed)1125 polygon = [[20,20], [40,20], [40,40], [20,40]]1126 points = [ [40, 50] ]1127 res = inside_polygon(points, polygon)1128 assert len(res) == 01129 1130 polygon = [[20,20], [40,20], [40,40], [20,40]]1131 points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]1132 res = inside_polygon(points, polygon)1133 assert len(res) == 21134 assert allclose(res, [0,1])1135 1136 1137 1138 #More convoluted and non convex polygon1139 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]1140 assert inside_polygon( (0.5, 0.5), polygon )1141 assert inside_polygon( (1, -0.5), polygon )1142 assert inside_polygon( (1.5, 0), polygon )1143 1144 assert not inside_polygon( (0.5, 1.5), polygon )1145 assert not inside_polygon( (0.5, -0.5), polygon )1146 1147 1148 #Very convoluted polygon1149 polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]1150 assert inside_polygon( (5, 5), polygon )1151 assert inside_polygon( (17, 7), polygon )1152 assert inside_polygon( (27, 2), polygon )1153 assert inside_polygon( (35, -5), polygon )1154 assert not inside_polygon( (15, 7), polygon )1155 assert not inside_polygon( (24, 3), polygon )1156 assert not inside_polygon( (25, -10), polygon )1157 1158 1159 1160 #Another combination (that failed)1161 polygon = [[0,0], [10,0], [10,10], [0,10]]1162 assert inside_polygon( (5, 5), polygon )1163 assert inside_polygon( (7, 7), polygon )1164 assert not inside_polygon( (-17, 7), polygon )1165 assert not inside_polygon( (7, 17), polygon )1166 assert not inside_polygon( (17, 7), polygon )1167 assert not inside_polygon( (27, 8), polygon )1168 assert not inside_polygon( (35, -5), polygon )1169 1170 1171 1172 1173 #Multiple polygons1174 1175 polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],1176 [10,10], [11,10], [11,11], [10,11], [10,10]]1177 assert inside_polygon( (0.5, 0.5), polygon )1178 assert inside_polygon( (10.5, 10.5), polygon )1179 1180 #FIXME: Fails if point is 5.5, 5.51181 assert not inside_polygon( (0, 5.5), polygon )1182 1183 #Polygon with a hole1184 polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],1185 [0,0], [1,0], [1,1], [0,1], [0,0]]1186 1187 assert inside_polygon( (0, -0.5), polygon )1188 assert not inside_polygon( (0.5, 0.5), polygon )1189 1190 def test_inside_polygon_vector_version(self):1191 #Now try the vector formulation returning indices1192 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]1193 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]1194 res = inside_polygon( points, polygon, verbose=False )1195 1196 assert allclose( res, [0,1,2] )1197 1198 def test_outside_polygon(self):1199 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square1200 1201 assert not outside_polygon( [0.5, 0.5], U )1202 #evaluate to False as the point 0.5, 0.5 is inside the unit square1203 1204 assert outside_polygon( [1.5, 0.5], U )1205 #evaluate to True as the point 1.5, 0.5 is outside the unit square1206 1207 indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )1208 assert allclose( indices, [1] )1209 1210 #One more test of vector formulation returning indices1211 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]1212 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]1213 res = outside_polygon( points, polygon )1214 1215 assert allclose( res, [3, 4] )1216 1217 1218 1219 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]1220 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]1221 res = outside_polygon( points, polygon )1222 1223 assert allclose( res, [0, 4, 5] )1224 1225 def test_outside_polygon2(self):1226 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square1227 1228 assert not outside_polygon( [0.5, 1.0], U, closed = True )1229 #evaluate to False as the point 0.5, 1.0 is inside the unit square1230 1231 assert outside_polygon( [0.5, 1.0], U, closed = False )1232 #evaluate to True as the point 0.5, 1.0 is outside the unit square1233 1234 def test_separate_points_by_polygon(self):1235 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square1236 1237 indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )1238 assert allclose( indices, [0,2,1] )1239 assert count == 21240 1241 #One more test of vector formulation returning indices1242 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]1243 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]1244 res, count = separate_points_by_polygon( points, polygon )1245 1246 assert allclose( res, [0,1,2,4,3] )1247 assert count == 31248 1249 1250 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]1251 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]1252 res, count = separate_points_by_polygon( points, polygon )1253 1254 assert allclose( res, [1,2,3,5,4,0] )1255 assert count == 31256 1257 1258 def test_populate_polygon(self):1259 1260 polygon = [[0,0], [1,0], [1,1], [0,1]]1261 points = populate_polygon(polygon, 5)1262 1263 assert len(points) == 51264 for point in points:1265 assert inside_polygon(point, polygon)1266 1267 1268 #Very convoluted polygon1269 polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]1270 1271 points = populate_polygon(polygon, 5)1272 1273 assert len(points) == 51274 for point in points:1275 assert inside_polygon(point, polygon)1276 1277 1278 def test_populate_polygon_with_exclude(self):1279 1280 1281 polygon = [[0,0], [1,0], [1,1], [0,1]]1282 ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter1283 points = populate_polygon(polygon, 5, exclude = [ex_poly])1284 1285 assert len(points) == 51286 for point in points:1287 assert inside_polygon(point, polygon)1288 assert not inside_polygon(point, ex_poly)1289 1290 1291 #overlap1292 polygon = [[0,0], [1,0], [1,1], [0,1]]1293 ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]1294 points = populate_polygon(polygon, 5, exclude = [ex_poly])1295 1296 assert len(points) == 51297 for point in points:1298 assert inside_polygon(point, polygon)1299 assert not inside_polygon(point, ex_poly)1300 1301 #Multiple1302 polygon = [[0,0], [1,0], [1,1], [0,1]]1303 ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter1304 ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter1305 1306 points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])1307 1308 assert len(points) == 201309 for point in points:1310 assert inside_polygon(point, polygon)1311 assert not inside_polygon(point, ex_poly1)1312 assert not inside_polygon(point, ex_poly2)1313 1314 1315 #Very convoluted polygon1316 polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]1317 ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]1318 points = populate_polygon(polygon, 20, exclude = [ex_poly])1319 1320 assert len(points) == 201321 for point in points:1322 assert inside_polygon(point, polygon)1323 assert not inside_polygon(point, ex_poly), '%s' %str(point)1324 975 1325 976 -
inundation/pyvolution/util.py
r1905 r1911 5 5 """ 6 6 7 #FIXME: Move polygon stuff out 7 #Soon to be obsoleted here 8 9 10 #FIXME: Obsolete this shortcut 11 from utilities.numerical_tools import ensure_numeric 12 13 14 import utilities.polygon 15 16 def point_on_line(*args, **kwargs): 17 """Temporary Interface to new location""" 18 19 print 'point_on_line has moved from util.py. ', 20 print 'Please use "from utilities.polygon import point_on_line"' 21 22 return utilities.polygon.point_on_line(*args, **kwargs) 23 24 def inside_polygon(*args, **kwargs): 25 """Temporary Interface to new location""" 26 27 print 'inside_polygon has moved from util.py. ', 28 print 'Please use "from utilities.polygon import inside_polygon"' 29 30 return utilities.polygon.inside_polygon(*args, **kwargs) 31 32 def outside_polygon(*args, **kwargs): 33 """Temporary Interface to new location""" 34 35 print 'outside_polygon has moved from util.py. ', 36 print 'Please use "from utilities.polygon import outside_polygon"' 37 38 return utilities.polygon.outside_polygon(*args, **kwargs) 39 40 41 def separate_points_by_polygon(*args, **kwargs): 42 """Temporary Interface to new location""" 43 44 print 'separate_points_by_polygon has moved from util.py. ', 45 print 'Please use "from utilities.polygon import separate_points_by_polygon"' 46 47 return utilities.polygon.separate_points_by_polygon(*args, **kwargs) 48 49 50 from utilities.polygon import Polygon_function #No warning 51 52 53 def read_polygon(*args, **kwargs): 54 """Temporary Interface to new location""" 55 56 print 'read_polygon has moved from util.py. ', 57 print 'Please use "from utilities.polygon import read_polygon"' 58 59 return utilities.polygon.read_polygon(*args, **kwargs) 60 61 62 def populate_polygon(*args, **kwargs): 63 """Temporary Interface to new location""" 64 65 print 'populate_polygon has moved from util.py. ', 66 print 'Please use "from utilities.polygon import populate_polygon"' 67 68 return utilities.polygon.populate_polygon(*args, **kwargs) 69 70 def read_xya(filename): 71 """Read simple xya file, possibly with a header in the first line, just 72 x y [attributes] 73 separated by whitespace 74 75 Format can be either ASCII or NetCDF 76 77 Return list of points, list of attributes and 78 attribute names if present otherwise None 79 """ 80 print "read_xya is obsolete. Use import_points_file from load_mesh.loadASCII" 81 #FIXME: Probably obsoleted by similar function in load_ASCII 82 #FIXME: Name read_data_points (or see 'load' in pylab) 83 84 85 from load_mesh.loadASCII import import_points_file 86 87 points_dict = import_points_file(filename) 88 return points_dict['pointlist'], points_dict['attributelist'] 89 90 91 92 93 #Normal util stuff 8 94 9 95 def angle(v): … … 73 159 74 160 75 def point_on_line(x, y, x0, y0, x1, y1):76 """Determine whether a point is on a line segment77 78 Input: x, y, x0, x0, x1, y1: where79 point is given by x, y80 line is given by (x0, y0) and (x1, y1)81 82 """83 84 from Numeric import array, dot, allclose85 from math import sqrt86 87 a = array([x - x0, y - y0])88 a_normal = array([a[1], -a[0]])89 90 b = array([x1 - x0, y1 - y0])91 92 if dot(a_normal, b) == 0:93 #Point is somewhere on the infinite extension of the line94 95 len_a = sqrt(sum(a**2))96 len_b = sqrt(sum(b**2))97 if dot(a, b) >= 0 and len_a <= len_b:98 return True99 else:100 return False101 else:102 return False103 104 def ensure_numeric(A, typecode = None):105 """Ensure that sequence is a Numeric array.106 Inputs:107 A: Sequance. If A is already a Numeric array it will be returned108 unaltered109 If not, an attempt is made to convert it to a Numeric110 array111 typecode: Numeric type. If specified, use this in the conversion.112 If not, let Numeric decide113 114 This function is necessary as array(A) can cause memory overflow.115 """116 117 from Numeric import ArrayType, array118 119 if typecode is None:120 if type(A) == ArrayType:121 return A122 else:123 return array(A)124 else:125 if type(A) == ArrayType:126 if A.typecode == typecode:127 return array(A) #FIXME: Shouldn't this be just return A?128 else:129 return A.astype(typecode)130 else:131 return array(A).astype(typecode)132 133 134 135 161 def file_function(filename, 136 162 domain = None, … … 393 419 394 420 395 def read_xya(filename): 396 """Read simple xya file, possibly with a header in the first line, just 397 x y [attributes] 398 separated by whitespace 399 400 Format can be either ASCII or NetCDF 401 402 Return list of points, list of attributes and 403 attribute names if present otherwise None 404 """ 405 print "read_xya is obsolete. Use import_points_file from load_mesh.loadASCII" 406 #FIXME: Probably obsoleted by similar function in load_ASCII 407 #FIXME: Name read_data_points (or see 'load' in pylab) 408 409 410 from load_mesh.loadASCII import import_points_file 411 412 points_dict = import_points_file(filename) 413 return points_dict['pointlist'], points_dict['attributelist'] 414 415 416 ##################################### 417 #POLYGON STUFF 418 # 419 #FIXME: All these should be put into new module polygon.py 420 421 422 423 #These have been redefined to use separate_by_polygon which will 424 #put all inside indices in the first part of the indices array and 425 #outside indices in the last 426 427 def inside_polygon(points, polygon, closed = True, verbose = False): 428 """See separate_points_by_polygon for documentation 429 """ 430 431 from Numeric import array, Float, reshape 432 433 if verbose: print 'Checking input to inside_polygon' 434 try: 435 points = ensure_numeric(points, Float) 436 except: 437 msg = 'Points could not be converted to Numeric array' 438 raise msg 439 440 try: 441 polygon = ensure_numeric(polygon, Float) 442 except: 443 msg = 'Polygon could not be converted to Numeric array' 444 raise msg 445 446 447 448 if len(points.shape) == 1: 449 one_point = True 450 points = reshape(points, (1,2)) 451 else: 452 one_point = False 453 454 indices, count = separate_points_by_polygon(points, polygon, 455 closed, verbose) 456 457 if one_point: 458 return count == 1 459 else: 460 return indices[:count] 461 462 def outside_polygon(points, polygon, closed = True, verbose = False): 463 """See separate_points_by_polygon for documentation 464 """ 465 466 from Numeric import array, Float, reshape 467 468 if verbose: print 'Checking input to outside_polygon' 469 try: 470 points = ensure_numeric(points, Float) 471 except: 472 msg = 'Points could not be converted to Numeric array' 473 raise msg 474 475 try: 476 polygon = ensure_numeric(polygon, Float) 477 except: 478 msg = 'Polygon could not be converted to Numeric array' 479 raise msg 480 481 482 483 if len(points.shape) == 1: 484 one_point = True 485 points = reshape(points, (1,2)) 486 else: 487 one_point = False 488 489 indices, count = separate_points_by_polygon(points, polygon, 490 closed, verbose) 491 492 if one_point: 493 return count != 1 494 else: 495 return indices[count:][::-1] #return reversed 496 497 498 def separate_points_by_polygon(points, polygon, 499 closed = True, verbose = False): 500 """Determine whether points are inside or outside a polygon 501 502 Input: 503 points - Tuple of (x, y) coordinates, or list of tuples 504 polygon - list of vertices of polygon 505 closed - (optional) determine whether points on boundary should be 506 regarded as belonging to the polygon (closed = True) 507 or not (closed = False) 508 509 Outputs: 510 indices: array of same length as points with indices of points falling 511 inside the polygon listed from the beginning and indices of points 512 falling outside listed from the end. 513 514 count: count of points falling inside the polygon 515 516 The indices of points inside are obtained as indices[:count] 517 The indices of points outside are obtained as indices[count:] 518 519 520 Examples: 521 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square 522 523 separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U) 524 will return the indices [0, 2, 1] and count == 2 as only the first 525 and the last point are inside the unit square 526 527 Remarks: 528 The vertices may be listed clockwise or counterclockwise and 529 the first point may optionally be repeated. 530 Polygons do not need to be convex. 531 Polygons can have holes in them and points inside a hole is 532 regarded as being outside the polygon. 533 534 Algorithm is based on work by Darel Finley, 535 http://www.alienryderflex.com/polygon/ 536 537 """ 538 539 from Numeric import array, Float, reshape, Int, zeros 540 541 542 #Input checks 543 try: 544 points = ensure_numeric(points, Float) 545 except: 546 msg = 'Points could not be converted to Numeric array' 547 raise msg 548 549 try: 550 polygon = ensure_numeric(polygon, Float) 551 except: 552 msg = 'Polygon could not be converted to Numeric array' 553 raise msg 554 555 assert len(polygon.shape) == 2,\ 556 'Polygon array must be a 2d array of vertices' 557 558 assert polygon.shape[1] == 2,\ 559 'Polygon array must have two columns' 560 561 assert len(points.shape) == 2,\ 562 'Points array must be a 2d array' 563 564 assert points.shape[1] == 2,\ 565 'Points array must have two columns' 566 567 N = polygon.shape[0] #Number of vertices in polygon 568 M = points.shape[0] #Number of points 569 570 px = polygon[:,0] 571 py = polygon[:,1] 572 573 #Used for an optimisation when points are far away from polygon 574 minpx = min(px); maxpx = max(px) 575 minpy = min(py); maxpy = max(py) 576 577 578 #Begin main loop 579 indices = zeros(M, Int) 580 581 inside_index = 0 #Keep track of points inside 582 outside_index = M-1 #Keep track of points outside (starting from end) 583 584 for k in range(M): 585 586 if verbose: 587 if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M) 588 589 #for each point 590 x = points[k, 0] 591 y = points[k, 1] 592 593 inside = False 594 595 if not x > maxpx or x < minpx or y > maxpy or y < minpy: 596 #Check polygon 597 for i in range(N): 598 j = (i+1)%N 599 600 #Check for case where point is contained in line segment 601 if point_on_line(x, y, px[i], py[i], px[j], py[j]): 602 if closed: 603 inside = True 604 else: 605 inside = False 606 break 607 else: 608 #Check if truly inside polygon 609 if py[i] < y and py[j] >= y or\ 610 py[j] < y and py[i] >= y: 611 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: 612 inside = not inside 613 614 if inside: 615 indices[inside_index] = k 616 inside_index += 1 617 else: 618 indices[outside_index] = k 619 outside_index -= 1 620 621 return indices, inside_index 622 623 624 def separate_points_by_polygon_c(points, polygon, 625 closed = True, verbose = False): 626 """Determine whether points are inside or outside a polygon 627 628 C-wrapper 629 """ 630 631 from Numeric import array, Float, reshape, zeros, Int 632 633 634 if verbose: print 'Checking input to separate_points_by_polygon' 635 #Input checks 636 try: 637 points = ensure_numeric(points, Float) 638 except: 639 msg = 'Points could not be converted to Numeric array' 640 raise msg 641 642 #if verbose: print 'Checking input to separate_points_by_polygon 2' 643 try: 644 polygon = ensure_numeric(polygon, Float) 645 except: 646 msg = 'Polygon could not be converted to Numeric array' 647 raise msg 648 649 if verbose: print 'check' 650 651 assert len(polygon.shape) == 2,\ 652 'Polygon array must be a 2d array of vertices' 653 654 assert polygon.shape[1] == 2,\ 655 'Polygon array must have two columns' 656 657 assert len(points.shape) == 2,\ 658 'Points array must be a 2d array' 659 660 assert points.shape[1] == 2,\ 661 'Points array must have two columns' 662 663 N = polygon.shape[0] #Number of vertices in polygon 664 M = points.shape[0] #Number of points 665 666 from util_ext import separate_points_by_polygon 667 668 if verbose: print 'Allocating array for indices' 669 670 indices = zeros( M, Int ) 671 672 if verbose: print 'Calling C-version of inside poly' 673 count = separate_points_by_polygon(points, polygon, indices, 674 int(closed), int(verbose)) 675 676 return indices, count 677 678 679 680 class Polygon_function: 681 """Create callable object f: x,y -> z, where a,y,z are vectors and 682 where f will return different values depending on whether x,y belongs 683 to specified polygons. 684 685 To instantiate: 686 687 Polygon_function(polygons) 688 689 where polygons is a list of tuples of the form 690 691 [ (P0, v0), (P1, v1), ...] 692 693 with Pi being lists of vertices defining polygons and vi either 694 constants or functions of x,y to be applied to points with the polygon. 695 696 The function takes an optional argument, default which is the value 697 (or function) to used for points not belonging to any polygon. 698 For example: 699 700 Polygon_function(polygons, default = 0.03) 701 702 If omitted the default value will be 0.0 703 704 Note: If two polygons overlap, the one last in the list takes precedence 705 706 FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference??? 707 708 """ 709 710 def __init__(self, regions, default = 0.0): 711 712 try: 713 len(regions) 714 except: 715 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 716 raise msg 717 718 719 T = regions[0] 720 try: 721 a = len(T) 722 except: 723 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 724 raise msg 725 726 assert a == 2, 'Must have two component each: %s' %T 727 728 self.regions = regions 729 self.default = default 730 731 732 def __call__(self, x, y): 733 from util import inside_polygon 734 from Numeric import ones, Float, concatenate, array, reshape, choose 735 736 x = array(x).astype(Float) 737 y = array(y).astype(Float) 738 739 N = len(x) 740 assert len(y) == N 741 742 points = concatenate( (reshape(x, (N, 1)), 743 reshape(y, (N, 1))), axis=1 ) 744 745 if callable(self.default): 746 z = self.default(x,y) 747 else: 748 z = ones(N, Float) * self.default 749 750 for polygon, value in self.regions: 751 indices = inside_polygon(points, polygon) 752 753 #FIXME: This needs to be vectorised 754 if callable(value): 755 for i in indices: 756 xx = array([x[i]]) 757 yy = array([y[i]]) 758 z[i] = value(xx, yy)[0] 759 else: 760 for i in indices: 761 z[i] = value 762 763 return z 764 765 def read_polygon(filename): 766 """Read points assumed to form a polygon 767 There must be exactly two numbers in each line 768 """ 769 770 #Get polygon 771 fid = open(filename) 772 lines = fid.readlines() 773 fid.close() 774 polygon = [] 775 for line in lines: 776 fields = line.split(',') 777 polygon.append( [float(fields[0]), float(fields[1])] ) 778 779 return polygon 780 781 def populate_polygon(polygon, number_of_points, seed = None, exclude = None): 782 """Populate given polygon with uniformly distributed points. 783 784 Input: 785 polygon - list of vertices of polygon 786 number_of_points - (optional) number of points 787 seed - seed for random number generator (default=None) 788 exclude - list of polygons (inside main polygon) from where points should be excluded 789 790 Output: 791 points - list of points inside polygon 792 793 Examples: 794 populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 ) 795 will return five randomly selected points inside the unit square 796 """ 797 798 from random import uniform, seed 799 800 seed(seed) 801 802 points = [] 803 804 #Find outer extent of polygon 805 max_x = min_x = polygon[0][0] 806 max_y = min_y = polygon[0][1] 807 for point in polygon[1:]: 808 x = point[0] 809 if x > max_x: max_x = x 810 if x < min_x: min_x = x 811 y = point[1] 812 if y > max_y: max_y = y 813 if y < min_y: min_y = y 814 815 816 while len(points) < number_of_points: 817 x = uniform(min_x, max_x) 818 y = uniform(min_y, max_y) 819 820 append = False 821 if inside_polygon( [x,y], polygon ): 822 823 append = True 824 825 #Check exclusions 826 if exclude is not None: 827 for ex_poly in exclude: 828 if inside_polygon( [x,y], ex_poly ): 829 append = False 830 831 832 if append is True: 833 points.append([x,y]) 834 835 return points 836 421 422 837 423 838 424 … … 914 500 import compile 915 501 if compile.can_use_C_extension('util_ext.c'): 916 from util_ext import gradient, gradient2 , point_on_line917 separate_points_by_polygon = separate_points_by_polygon_c502 from util_ext import gradient, gradient2#, point_on_line 503 #separate_points_by_polygon = separate_points_by_polygon_c 918 504 else: 919 505 gradient = gradient_python -
inundation/utilities/polygon.py
r1910 r1911 36 36 37 37 38 39 #These have been redefined to use separate_by_polygon which will40 #put all inside indices in the first part of the indices array and41 #outside indices in the last42 43 38 def inside_polygon(points, polygon, closed = True, verbose = False): 44 """See separate_points_by_polygon for documentation 39 """Determine points inside a polygon 40 41 Functions inside_polygon and outside_polygon have been defined in 42 terms af separate_by_polygon which will put all inside indices in 43 the first part of the indices array and outside indices in the last 44 45 See separate_points_by_polygon for documentation 45 46 """ 46 47 … … 77 78 78 79 def outside_polygon(points, polygon, closed = True, verbose = False): 79 """See separate_points_by_polygon for documentation 80 """Determine points outside a polygon 81 82 Functions inside_polygon and outside_polygon have been defined in 83 terms af separate_by_polygon which will put all inside indices in 84 the first part of the indices array and outside indices in the last 85 86 See separate_points_by_polygon for documentation 80 87 """ 81 88
Note: See TracChangeset
for help on using the changeset viewer.