source:inundation/utilities/polygon.py@1906

Last change on this file since 1906 was 1906, checked in by ole, 18 years ago

Began moving polygon functionality to separate module and created a new package utilities for general tools used by AnuGa?

File size: 12.7 KB
Line
1
2
3
4
5#####################################
6#POLYGON STUFF
7#
8#FIXME: All these should be put into new module polygon.py
9
10
11
12def point_on_line(x, y, x0, y0, x1, y1):
13    """Determine whether a point is on a line segment
14
15    Input: x, y, x0, x0, x1, y1: where
16        point is given by x, y
17        line is given by (x0, y0) and (x1, y1)
18
19    """
20
21    from Numeric import array, dot, allclose
22    from math import sqrt
23
24    a = array([x - x0, y - y0])
25    a_normal = array([a[1], -a[0]])
26
27    b = array([x1 - x0, y1 - y0])
28
29    if dot(a_normal, b) == 0:
30        #Point is somewhere on the infinite extension of the line
31
32        len_a = sqrt(sum(a**2))
33        len_b = sqrt(sum(b**2))
34        if dot(a, b) >= 0 and len_a <= len_b:
35           return True
36        else:
37           return False
38    else:
39      return False
40
41
42
43#These have been redefined to use separate_by_polygon which will
44#put all inside indices in the first part of the indices array and
45#outside indices in the last
46
47def inside_polygon(points, polygon, closed = True, verbose = False):
48    """See separate_points_by_polygon for documentation
49    """
50
51    from Numeric import array, Float, reshape
52
53    if verbose: print 'Checking input to inside_polygon'
54    try:
55        points = ensure_numeric(points, Float)
56    except:
57        msg = 'Points could not be converted to Numeric array'
58        raise msg
59
60    try:
61        polygon = ensure_numeric(polygon, Float)
62    except:
63        msg = 'Polygon could not be converted to Numeric array'
64        raise msg
65
66
67
68    if len(points.shape) == 1:
69        one_point = True
70        points = reshape(points, (1,2))
71    else:
72        one_point = False
73
74    indices, count = separate_points_by_polygon(points, polygon,
75                                                closed, verbose)
76
77    if one_point:
78        return count == 1
79    else:
80        return indices[:count]
81
82def outside_polygon(points, polygon, closed = True, verbose = False):
83    """See separate_points_by_polygon for documentation
84    """
85
86    from Numeric import array, Float, reshape
87
88    if verbose: print 'Checking input to outside_polygon'
89    try:
90        points = ensure_numeric(points, Float)
91    except:
92        msg = 'Points could not be converted to Numeric array'
93        raise msg
94
95    try:
96        polygon = ensure_numeric(polygon, Float)
97    except:
98        msg = 'Polygon could not be converted to Numeric array'
99        raise msg
100
101
102
103    if len(points.shape) == 1:
104        one_point = True
105        points = reshape(points, (1,2))
106    else:
107        one_point = False
108
109    indices, count = separate_points_by_polygon(points, polygon,
110                                                closed, verbose)
111
112    if one_point:
113        return count != 1
114    else:
115        return indices[count:][::-1]  #return reversed
116
117
118def separate_points_by_polygon(points, polygon,
119                               closed = True, verbose = False):
120    """Determine whether points are inside or outside a polygon
121
122    Input:
123       points - Tuple of (x, y) coordinates, or list of tuples
124       polygon - list of vertices of polygon
125       closed - (optional) determine whether points on boundary should be
126       regarded as belonging to the polygon (closed = True)
127       or not (closed = False)
128
129    Outputs:
130       indices: array of same length as points with indices of points falling
131       inside the polygon listed from the beginning and indices of points
132       falling outside listed from the end.
133
134       count: count of points falling inside the polygon
135
136       The indices of points inside are obtained as indices[:count]
137       The indices of points outside are obtained as indices[count:]
138
139
140    Examples:
141       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
142
143       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
144       will return the indices [0, 2, 1] and count == 2 as only the first
145       and the last point are inside the unit square
146
147    Remarks:
148       The vertices may be listed clockwise or counterclockwise and
149       the first point may optionally be repeated.
150       Polygons do not need to be convex.
151       Polygons can have holes in them and points inside a hole is
152       regarded as being outside the polygon.
153
154    Algorithm is based on work by Darel Finley,
155    http://www.alienryderflex.com/polygon/
156
157    """
158
159    from Numeric import array, Float, reshape, Int, zeros
160
161
162    #Input checks
163    try:
164        points = ensure_numeric(points, Float)
165    except:
166        msg = 'Points could not be converted to Numeric array'
167        raise msg
168
169    try:
170        polygon = ensure_numeric(polygon, Float)
171    except:
172        msg = 'Polygon could not be converted to Numeric array'
173        raise msg
174
175    assert len(polygon.shape) == 2,\
176       'Polygon array must be a 2d array of vertices'
177
178    assert polygon.shape[1] == 2,\
179       'Polygon array must have two columns'
180
181    assert len(points.shape) == 2,\
182       'Points array must be a 2d array'
183
184    assert points.shape[1] == 2,\
185       'Points array must have two columns'
186
187    N = polygon.shape[0] #Number of vertices in polygon
188    M = points.shape[0]  #Number of points
189
190    px = polygon[:,0]
191    py = polygon[:,1]
192
193    #Used for an optimisation when points are far away from polygon
194    minpx = min(px); maxpx = max(px)
195    minpy = min(py); maxpy = max(py)
196
197
198    #Begin main loop
199    indices = zeros(M, Int)
200
201    inside_index = 0    #Keep track of points inside
202    outside_index = M-1 #Keep track of points outside (starting from end)
203
204    for k in range(M):
205
206        if verbose:
207            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
208
209        #for each point
210        x = points[k, 0]
211        y = points[k, 1]
212
213        inside = False
214
215        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
216            #Check polygon
217            for i in range(N):
218                j = (i+1)%N
219
220                #Check for case where point is contained in line segment
221                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
222                    if closed:
223                        inside = True
224                    else:
225                        inside = False
226                    break
227                else:
228                    #Check if truly inside polygon
229                    if py[i] < y and py[j] >= y or\
230                       py[j] < y and py[i] >= y:
231                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
232                            inside = not inside
233
234        if inside:
235            indices[inside_index] = k
236            inside_index += 1
237        else:
238            indices[outside_index] = k
239            outside_index -= 1
240
241    return indices, inside_index
242
243
244def separate_points_by_polygon_c(points, polygon,
245                                 closed = True, verbose = False):
246    """Determine whether points are inside or outside a polygon
247
248    C-wrapper
249    """
250
251    from Numeric import array, Float, reshape, zeros, Int
252
253
254    if verbose: print 'Checking input to separate_points_by_polygon'
255    #Input checks
256    try:
257        points = ensure_numeric(points, Float)
258    except:
259        msg = 'Points could not be converted to Numeric array'
260        raise msg
261
262    #if verbose: print 'Checking input to separate_points_by_polygon 2'
263    try:
264        polygon = ensure_numeric(polygon, Float)
265    except:
266        msg = 'Polygon could not be converted to Numeric array'
267        raise msg
268
269    if verbose: print 'check'
270
271    assert len(polygon.shape) == 2,\
272       'Polygon array must be a 2d array of vertices'
273
274    assert polygon.shape[1] == 2,\
275       'Polygon array must have two columns'
276
277    assert len(points.shape) == 2,\
278       'Points array must be a 2d array'
279
280    assert points.shape[1] == 2,\
281       'Points array must have two columns'
282
283    N = polygon.shape[0] #Number of vertices in polygon
284    M = points.shape[0]  #Number of points
285
286    from util_ext import separate_points_by_polygon
287
288    if verbose: print 'Allocating array for indices'
289
290    indices = zeros( M, Int )
291
292    if verbose: print 'Calling C-version of inside poly'
293    count = separate_points_by_polygon(points, polygon, indices,
294                                       int(closed), int(verbose))
295
296    return indices, count
297
298
299
300class Polygon_function:
301    """Create callable object f: x,y -> z, where a,y,z are vectors and
302    where f will return different values depending on whether x,y belongs
303    to specified polygons.
304
305    To instantiate:
306
307       Polygon_function(polygons)
308
309    where polygons is a list of tuples of the form
310
311      [ (P0, v0), (P1, v1), ...]
312
313      with Pi being lists of vertices defining polygons and vi either
314      constants or functions of x,y to be applied to points with the polygon.
315
316    The function takes an optional argument, default which is the value
317    (or function) to used for points not belonging to any polygon.
318    For example:
319
320       Polygon_function(polygons, default = 0.03)
321
322    If omitted the default value will be 0.0
323
324    Note: If two polygons overlap, the one last in the list takes precedence
325
326    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???
327
328    """
329
330    def __init__(self, regions, default = 0.0):
331
332        try:
333            len(regions)
334        except:
335            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
336            raise msg
337
338
339        T = regions[0]
340        try:
341            a = len(T)
342        except:
343            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
344            raise msg
345
346        assert a == 2, 'Must have two component each: %s' %T
347
348        self.regions = regions
349        self.default = default
350
351
352    def __call__(self, x, y):
353        from util import inside_polygon
354        from Numeric import ones, Float, concatenate, array, reshape, choose
355
356        x = array(x).astype(Float)
357        y = array(y).astype(Float)
358
359        N = len(x)
360        assert len(y) == N
361
362        points = concatenate( (reshape(x, (N, 1)),
363                               reshape(y, (N, 1))), axis=1 )
364
365        if callable(self.default):
366            z = self.default(x,y)
367        else:
368            z = ones(N, Float) * self.default
369
370        for polygon, value in self.regions:
371            indices = inside_polygon(points, polygon)
372
373            #FIXME: This needs to be vectorised
374            if callable(value):
375                for i in indices:
376                    xx = array([x[i]])
377                    yy = array([y[i]])
378                    z[i] = value(xx, yy)[0]
379            else:
380                for i in indices:
381                    z[i] = value
382
383        return z
384
386    """Read points assumed to form a polygon
387       There must be exactly two numbers in each line
388    """
389
390    #Get polygon
391    fid = open(filename)
393    fid.close()
394    polygon = []
395    for line in lines:
396        fields = line.split(',')
397        polygon.append( [float(fields[0]), float(fields[1])] )
398
399    return polygon
400
401def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
402    """Populate given polygon with uniformly distributed points.
403
404    Input:
405       polygon - list of vertices of polygon
406       number_of_points - (optional) number of points
407       seed - seed for random number generator (default=None)
408       exclude - list of polygons (inside main polygon) from where points should be excluded
409
410    Output:
411       points - list of points inside polygon
412
413    Examples:
414       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
415       will return five randomly selected points inside the unit square
416    """
417
418    from random import uniform, seed
419
420    seed(seed)
421
422    points = []
423
424    #Find outer extent of polygon
425    max_x = min_x = polygon[0][0]
426    max_y = min_y = polygon[0][1]
427    for point in polygon[1:]:
428        x = point[0]
429        if x > max_x: max_x = x
430        if x < min_x: min_x = x
431        y = point[1]
432        if y > max_y: max_y = y
433        if y < min_y: min_y = y
434
435
436    while len(points) < number_of_points:
437        x = uniform(min_x, max_x)
438        y = uniform(min_y, max_y)
439
440        append = False
441        if inside_polygon( [x,y], polygon ):
442
443            append = True
444
445            #Check exclusions
446            if exclude is not None:
447                for ex_poly in exclude:
448                    if inside_polygon( [x,y], ex_poly ):
449                        append = False
450
451
452        if append is True:
453            points.append([x,y])
454
455    return points
456
457
458
459
460
461##############################################
462#Initialise module
463
464import compile
465if compile.can_use_C_extension('polygon_ext.c'):
466    from polygon_ext import point_on_line
467    separate_points_by_polygon = separate_points_by_polygon_c
468
469
470if __name__ == "__main__":
471    pass
Note: See TracBrowser for help on using the repository browser.