source: inundation/utilities/polygon.py @ 1911

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

Refactored pyvolution to use polygon functionality from new utilities package

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