source: trunk/anuga_core/source/anuga/geometry/polygon_function.py @ 7876

Last change on this file since 7876 was 7858, checked in by hudson, 14 years ago

Refactorings to increase code quality, fixed missing log import from sww_interrogate.

File size: 4.7 KB
Line 
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
7import anuga.utilities.log as log
8import numpy as num
9from polygon import inside_polygon
10
11class 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
Note: See TracBrowser for help on using the repository browser.