Changeset 7778 for trunk/anuga_core/source/anuga/geometry/polygon.py
- Timestamp:
- Jun 4, 2010, 10:03:48 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/geometry/polygon.py
r7711 r7778 7 7 from math import sqrt 8 8 from anuga.utilities.numerical_tools import ensure_numeric 9 from anuga.geospatial_data.geospatial_data import ensure_absolute, Geospatial_data 9 from anuga.geospatial_data.geospatial_data import ensure_absolute, \ 10 Geospatial_data 10 11 from anuga.config import netcdf_float 11 12 import anuga.utilities.log as log … … 69 70 % (str(p1), str(p2), str(p3), str(p4))) 70 71 71 # 0s1 0e1 1s0 1e0 # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0 72 collinear_result = { (False, False, False, False): lines_dont_coincide, 73 (False, False, False, True ): lines_error, 74 (False, False, True, False): lines_error, 75 (False, False, True, True ): lines_1_fully_included_in_0, 76 (False, True, False, False): lines_error, 77 (False, True, False, True ): lines_overlap_opposite_direction2, 78 (False, True, True, False): lines_overlap_same_direction2, 79 (False, True, True, True ): lines_1_fully_included_in_0, 80 (True, False, False, False): lines_error, 81 (True, False, False, True ): lines_overlap_same_direction, 82 (True, False, True, False): lines_overlap_opposite_direction, 83 (True, False, True, True ): lines_1_fully_included_in_0, 84 (True, True, False, False): lines_0_fully_included_in_1, 85 (True, True, False, True ): lines_0_fully_included_in_1, 86 (True, True, True, False): lines_0_fully_included_in_1, 87 (True, True, True, True ): lines_0_fully_included_in_1 88 } 72 collinear_result = { 73 # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0 74 # 0s1 0e1 1s0 1e0 75 (False, False, False, False): lines_dont_coincide, 76 (False, False, False, True ): lines_error, 77 (False, False, True, False): lines_error, 78 (False, False, True, True ): lines_1_fully_included_in_0, 79 (False, True, False, False): lines_error, 80 (False, True, False, True ): lines_overlap_opposite_direction2, 81 (False, True, True, False): lines_overlap_same_direction2, 82 (False, True, True, True ): lines_1_fully_included_in_0, 83 (True, False, False, False): lines_error, 84 (True, False, False, True ): lines_overlap_same_direction, 85 (True, False, True, False): lines_overlap_opposite_direction, 86 (True, False, True, True ): lines_1_fully_included_in_0, 87 (True, True, False, False): lines_0_fully_included_in_1, 88 (True, True, False, True ): lines_0_fully_included_in_1, 89 (True, True, True, False): lines_0_fully_included_in_1, 90 (True, True, True, True ): lines_0_fully_included_in_1 91 } 89 92 90 93 ## … … 275 278 276 279 # Quickly reject points that are clearly outside 277 if point[0] < min(triangle[:, 0]): return False278 if point[0] > max(triangle[:, 0]): return False279 if point[1] < min(triangle[:, 1]): return False280 if point[1] > max(triangle[:, 1]): return False280 if point[0] < min(triangle[:, 0]): return False 281 if point[0] > max(triangle[:, 0]): return False 282 if point[1] < min(triangle[:, 1]): return False 283 if point[1] > max(triangle[:, 1]): return False 281 284 282 285 … … 367 370 (type, point) = intersection(leftmost, l_x[cmp]) 368 371 comparisons += 1 369 if type != 0 and type != 4 or (type == 2 and list(point[0]) != list(point[1])): 372 if type != 0 and type != 4 or (type == 2 and list(point[0]) !=\ 373 list(point[1])): 370 374 if verbose: 371 print 'Self-intersecting polygon found, type ', type, ' point', point, 375 print 'Self-intersecting polygon found, type ', type 376 print 'point', point, 372 377 print 'vertices: ', leftmost, ' - ', l_x[cmp] 373 378 return True … … 586 591 if count == len(indices): 587 592 # No points are outside 588 return indices[:count], []593 return indices[:count], [] 589 594 else: 590 595 return indices[:count], indices[count:][::-1] #return reversed 591 596 592 ## 593 # @brief Sort a list of points into contiguous points inside+outside a polygon. 594 # @param points A set of points (tuple, list or array). 595 # @param polygon A set of points defining a polygon (tuple, list or array). 596 # @param closed True if points on boundary are considered 'inside' polygon. 597 # @param verbose True if this function is to be verbose. 598 # @return (indices, count) where indices are point indices and count is the 599 # delimiter index between point inside (on left) and others. 597 598 600 599 def separate_points_by_polygon(points, polygon, 601 600 closed=True, … … 644 643 if check_input: 645 644 #Input checks 646 assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean' 647 assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean' 645 assert isinstance(closed, bool), \ 646 'Keyword argument "closed" must be boolean' 647 assert isinstance(verbose, bool), \ 648 'Keyword argument "verbose" must be boolean' 648 649 649 650 try: … … 737 738 return abs(poly_area/2) 738 739 739 ## 740 # @brief Plot a set of polygons. 741 # @param polygons_points List of polygons to plot. 742 # @param style List of styles for each polygon. 743 # @param figname Name to save figure to. 744 # @param label Title for the plot. 745 # @param verbose True if this function is to be verbose. 746 # @return A list of min/max x and y values [minx, maxx, miny, maxy]. 747 # @note A style value is 'line' for polygons, 'outside' for points outside. 740 748 741 def plot_polygons(polygons_points, 749 742 style=None, … … 818 811 if s == 'outside': colour.append('r.') 819 812 if s == 'point': colour.append('g.') 820 if s <>'line':821 if s <>'outside':822 if s <>'point':813 if s != 'line': 814 if s != 'outside': 815 if s != 'point': 823 816 colour.append(s) 824 817 … … 836 829 title(label) 837 830 838 #raw_input('wait 1')839 #FIXME(Ole): This makes for some strange scalings sometimes.840 #if minx <> 0:841 # axis([minx*0.9,maxx*1.1,miny*0.9,maxy*1.1])842 #else:843 # if miny == 0:844 # axis([-maxx*.01,maxx*1.1,-maxy*0.01,maxy*1.1])845 # else:846 # axis([-maxx*.01,maxx*1.1,miny*0.9,maxy*1.1])847 848 831 if figname is not None: 849 832 savefig(figname) … … 856 839 return vec 857 840 858 ## 859 # @brief 860 # @param polygon A set of points defining a polygon. 861 # @param verbose True if this function is to be verbose. 862 # @return A tuple (x, y) of X and Y coordinates of the polygon. 863 # @note We duplicate the first point so can have closed polygon in plot. 841 864 842 def poly_xy(polygon, verbose=False): 865 843 """ this is used within plot_polygons so need to duplicate 866 844 the first point so can have closed polygon in plot 845 # @param polygon A set of points defining a polygon. 846 # @param verbose True if this function is to be verbose. 847 # @return A tuple (x, y) of X and Y coordinates of the polygon. 848 # @note We duplicate the first point so can have closed polygon in plot. 867 849 """ 868 850 … … 884 866 885 867 886 ##887 # @brief Define a class that defines a callable object for a polygon.888 # @note Object created is function: f: x,y -> z889 # where x, y and z are vectors and z depends on whether x,y belongs890 # to specified polygons.891 868 class Polygon_function: 892 869 """Create callable object f: x,y -> z, where a,y,z are vectors and … … 923 900 """ 924 901 925 ##926 # @brief Create instance of a polygon function.927 # @param regions A list of (x,y) tuples defining a polygon.928 # @param default Value or function returning value for points outside poly.929 # @param geo_reference ??930 902 def __init__(self, regions, default=0.0, geo_reference=None): 903 """ 904 # @brief Create instance of a polygon function. 905 # @param regions A list of (x,y) tuples defining a polygon. 906 # @param default Value or function returning value for points outside poly. 907 # @param geo_reference ?? 908 """ 909 931 910 try: 932 911 len(regions) … … 967 946 self.regions.append((P, value)) 968 947 969 ##970 # @brief Implement the 'callable' property of Polygon_function.971 # @param x List of x coordinates of points ot interest.972 # @param y List of y coordinates of points ot interest.973 948 def __call__(self, x, y): 949 """ 950 # @brief Implement the 'callable' property of Polygon_function. 951 # @param x List of x coordinates of points ot interest. 952 # @param y List of y coordinates of points ot interest. 953 """ 974 954 x = num.array(x, num.float) 975 955 y = num.array(y, num.float) … … 1014 994 ################################################################################ 1015 995 1016 ## 1017 # @brief Read polygon data from a file. 996 def read_polygon(filename, delimiter=','): 997 """Read points assumed to form a polygon. 998 1018 999 # @param filename Path to file containing polygon data. 1019 1000 # @param delimiter Delimiter to split polygon data with. 1020 1001 # @return A list of point data from the polygon file. 1021 def read_polygon(filename, delimiter=','): 1022 """Read points assumed to form a polygon. 1002 1023 1003 1024 1004 There must be exactly two numbers in each line separated by the delimiter. … … 1036 1016 # check this is a valid polygon. 1037 1017 if is_complex(polygon, verbose=True): 1038 msg = 'ERROR: Self-intersecting polygon detected in file ' + filename +'. ' 1039 msg += 'A complex polygon will not necessarily break the algorithms within ANUGA, ' 1040 msg += 'but it usually signifies pathological data. Please fix this file.' 1018 msg = 'ERROR: Self-intersecting polygon detected in file ' 1019 msg += filename +'. A complex polygon will not ' 1020 msg += 'necessarily break the algorithms within ANUGA, but it' 1021 msg += 'usually signifies pathological data. Please fix this file.' 1041 1022 raise Exception, msg 1042 1023 … … 1069 1050 pass 1070 1051 1071 ## 1072 # @brief Populate given polygon with uniformly distributed points. 1073 # @param polygon Polygon to uniformly fill. 1074 # @param number_of_points Number of points required in polygon. 1075 # @param seed Seed for random number generator. 1076 # @param exclude List of polygons inside main where points should be excluded. 1077 # @return List of random points inside input polygon. 1078 # @note Delimiter is assumed to be a comma. 1052 1079 1053 def populate_polygon(polygon, number_of_points, seed=None, exclude=None): 1080 1054 """Populate given polygon with uniformly distributed points. … … 1131 1105 return points 1132 1106 1133 ## 1134 # @brief Get a point inside a polygon that is close to an edge. 1135 # @param polygon List of vertices of polygon. 1136 # @param delta Maximum distance from an edge is delta * sqrt(2). 1137 # @return A point that is inside polgon and close to the polygon edge. 1107 1138 1108 def point_in_polygon(polygon, delta=1e-8): 1139 1109 """Return a point inside a given polygon which will be close to the … … 1186 1156 return point 1187 1157 1188 ## 1189 # @brief Calculate approximate number of triangles inside a bounding polygon. 1190 # @param interior_regions 1191 # @param bounding_poly 1192 # @param remainder_res 1193 # @return The number of triangles. 1158 1194 1159 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res): 1195 1160 """Calculate the approximate number of triangles inside the … … 1247 1212 return int(total_number_of_triangles) 1248 1213 1214 1215 def decimate_polygon(polygon, factor=10): 1216 """Reduce number of points in polygon by the specified 1217 factor (default=10, hence the name of the function) such that 1218 the extrema in both axes are preserved. 1219 1249 1220 ## 1250 1221 # @brief Reduce number of points in polygon by the specified factor. 1251 1222 # @param polygon The polygon to reduce. 1252 1223 # @param factor The factor to reduce polygon points by (default 10). 1253 # @return The reduced polygon points list.1254 1224 # @note The extrema of both axes are preserved. 1255 def decimate_polygon(polygon, factor=10):1256 """Reduce number of points in polygon by the specified1257 factor (default=10, hence the name of the function) such that1258 the extrema in both axes are preserved.1259 1225 1260 1226 Return reduced polygon … … 1285 1251 return reduced_polygon 1286 1252 1287 ## 1288 # @brief Interpolate linearly from polyline nodes to midpoints of triangles. 1289 # @param data The data on the polyline nodes. 1290 # @param polyline_nodes ?? 1291 # @param gauge_neighbour_id ?? FIXME(Ole): I want to get rid of this 1292 # @param point_coordinates ?? 1293 # @param verbose True if this function is to be verbose. 1253 1294 1254 def interpolate_polyline(data, 1295 1255 polyline_nodes, … … 1370 1330 1371 1331 return points, vertices 1372 ##1373 # @brief1374 # @param data1375 # @param polyline_nodes1376 # @param gauge_neighbour_id1377 # @param interpolation_points1378 # @param interpolated_values1379 # @param rtol1380 # @param atol1381 # @return1382 # @note OBSOLETED BY C-EXTENSION1383 def _interpolate_polyline(data,1384 polyline_nodes,1385 gauge_neighbour_id,1386 interpolation_points,1387 interpolated_values,1388 rtol=1.0e-6,1389 atol=1.0e-8):1390 """Auxiliary function used by interpolate_polyline1391 1392 NOTE: OBSOLETED BY C-EXTENSION1393 """1394 1395 number_of_nodes = len(polyline_nodes)1396 number_of_points = len(interpolation_points)1397 1398 for j in range(number_of_nodes):1399 neighbour_id = gauge_neighbour_id[j]1400 1401 # FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded,1402 # but need to check with John J.1403 # Keep it for now (17 Jan 2009)1404 # When gone, we can simply interpolate between neighbouring nodes,1405 # i.e. neighbour_id = j+1.1406 # and the test below becomes something like: if j < number_of_nodes...1407 1408 if neighbour_id >= 0:1409 x0, y0 = polyline_nodes[j,:]1410 x1, y1 = polyline_nodes[neighbour_id,:]1411 1412 segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)1413 segment_delta = data[neighbour_id] - data[j]1414 slope = segment_delta/segment_len1415 1416 for i in range(number_of_points):1417 x, y = interpolation_points[i,:]1418 if point_on_line([x, y], [[x0, y0], [x1, y1]],1419 rtol=rtol, atol=atol):1420 dist = sqrt((x-x0)**2 + (y-y0)**2)1421 interpolated_values[i] = slope*dist + data[j]1422 1332 1423 1333
Note: See TracChangeset
for help on using the changeset viewer.
