1 | """ |
---|
2 | Callable function to determine if points lie inside or outside a polygon. |
---|
3 | |
---|
4 | As of June 2010 this module has a pylint quality rating of 8.85/10. |
---|
5 | """ |
---|
6 | |
---|
7 | import anuga.utilities.log as log |
---|
8 | import numpy as num |
---|
9 | from polygon import inside_polygon |
---|
10 | |
---|
11 | class Polygon_function: |
---|
12 | """Create callable object f: x,y -> z, where a,y,z are vectors and |
---|
13 | where f will return different values depending on whether x,y belongs |
---|
14 | to specified polygons. |
---|
15 | |
---|
16 | To instantiate: |
---|
17 | |
---|
18 | Polygon_function(polygons) |
---|
19 | |
---|
20 | where polygons is a list of tuples of the form |
---|
21 | |
---|
22 | [ (P0, v0), (P1, v1), ...] |
---|
23 | |
---|
24 | with Pi being lists of vertices defining polygons and vi either |
---|
25 | constants or functions of x,y to be applied to points with the polygon. |
---|
26 | |
---|
27 | The function takes an optional argument, default which is the value |
---|
28 | (or function) to used for points not belonging to any polygon. |
---|
29 | For example: |
---|
30 | |
---|
31 | Polygon_function(polygons, default = 0.03) |
---|
32 | |
---|
33 | If omitted the default value will be 0.0 |
---|
34 | |
---|
35 | Note: If two polygons overlap, the one last in the list takes precedence |
---|
36 | |
---|
37 | Coordinates specified in the call are assumed to be relative to the |
---|
38 | origin (georeference) e.g. used by domain. |
---|
39 | By specifying the optional argument georeference, |
---|
40 | all points are made relative. |
---|
41 | |
---|
42 | FIXME: This should really work with geo_spatial point sets. |
---|
43 | """ |
---|
44 | |
---|
45 | def __init__(self, regions, default=0.0, geo_reference=None): |
---|
46 | """ |
---|
47 | # @brief Create instance of a polygon function. |
---|
48 | # @param regions A list of (x,y) tuples defining a polygon. |
---|
49 | # @param default Value or function returning value for points outside poly. |
---|
50 | # @param geo_reference ?? |
---|
51 | """ |
---|
52 | |
---|
53 | try: |
---|
54 | len(regions) |
---|
55 | except: |
---|
56 | msg = ('Polygon_function takes a list of pairs (polygon, value).' |
---|
57 | 'Got %s' % str(regions)) |
---|
58 | raise Exception, msg |
---|
59 | |
---|
60 | first_region = regions[0] |
---|
61 | |
---|
62 | if isinstance(first_region, basestring): |
---|
63 | msg = ('You passed in a list of text values into polygon_function ' |
---|
64 | 'instead of a list of pairs (polygon, value): "%s"' |
---|
65 | % str(first_region)) |
---|
66 | raise Exception, msg |
---|
67 | |
---|
68 | try: |
---|
69 | num_region_components = len(first_region) |
---|
70 | except: |
---|
71 | msg = ('Polygon_function takes a list of pairs (polygon, value). ' |
---|
72 | 'Got %s' % str(num_region_components)) |
---|
73 | raise Exception, msg |
---|
74 | |
---|
75 | msg = ('Each entry in regions have two components: (polygon, value). ' |
---|
76 | 'I got %s' % str(num_region_components)) |
---|
77 | assert num_region_components == 2, msg |
---|
78 | |
---|
79 | if geo_reference is None: |
---|
80 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
81 | geo_reference = Geo_reference() |
---|
82 | |
---|
83 | self.default = default |
---|
84 | |
---|
85 | # Make points in polygons relative to geo_reference |
---|
86 | self.regions = [] |
---|
87 | for polygon, value in regions: |
---|
88 | georeffed_poly = geo_reference.change_points_geo_ref(polygon) |
---|
89 | self.regions.append((georeffed_poly, value)) |
---|
90 | |
---|
91 | def __call__(self, pts_x, pts_y): |
---|
92 | """ |
---|
93 | # @brief Implement the 'callable' property of Polygon_function. |
---|
94 | # @param x List of x coordinates of points ot interest. |
---|
95 | # @param y List of y coordinates of points ot interest. |
---|
96 | """ |
---|
97 | pts_x = num.array(pts_x, num.float) |
---|
98 | pts_y = num.array(pts_y, num.float) |
---|
99 | |
---|
100 | # x and y must be one-dimensional and same length |
---|
101 | assert len(pts_x.shape) == 1 and len(pts_y.shape) == 1 |
---|
102 | pts_len = pts_x.shape[0] |
---|
103 | assert pts_y.shape[0] == pts_len, 'x and y must be same length' |
---|
104 | |
---|
105 | points = num.ascontiguousarray(num.concatenate((pts_x[:, num.newaxis], |
---|
106 | pts_y[:, num.newaxis]), |
---|
107 | axis = 1 )) |
---|
108 | |
---|
109 | if callable(self.default): |
---|
110 | result = self.default(pts_x, pts_y) |
---|
111 | else: |
---|
112 | result = num.ones(pts_len, num.float) * self.default |
---|
113 | |
---|
114 | for polygon, value in self.regions: |
---|
115 | indices = inside_polygon(points, polygon) |
---|
116 | |
---|
117 | # FIXME: This needs to be vectorised |
---|
118 | if callable(value): |
---|
119 | for i in indices: |
---|
120 | xx = num.array([pts_x[i]]) |
---|
121 | yy = num.array([pts_y[i]]) |
---|
122 | result[i] = value(xx, yy)[0] |
---|
123 | else: |
---|
124 | for i in indices: |
---|
125 | result[i] = value |
---|
126 | |
---|
127 | if len(result) == 0: |
---|
128 | msg = ('Warning: points provided to Polygon function did not fall ' |
---|
129 | 'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]' |
---|
130 | % (min(pts_x), max(pts_x), min(pts_y), max(pts_y))) |
---|
131 | log.critical(msg) |
---|
132 | |
---|
133 | return result |
---|