source:inundation/utilities/polygon.py@2309

Last change on this file since 2309 was 2241, checked in by duncan, 17 years ago

bug fix

File size: 14.2 KB
Line
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):
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
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):
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
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        if count == len(indices):
119            # No points are outside
120            return []
121        else:
122            return indices[count:][::-1]  #return reversed
123
124
125def separate_points_by_polygon(points, polygon,
126                               closed = True, verbose = False):
127    """Determine whether points are inside or outside a polygon
128
129    Input:
130       points - Tuple of (x, y) coordinates, or list of tuples
131       polygon - list of vertices of polygon
132       closed - (optional) determine whether points on boundary should be
133       regarded as belonging to the polygon (closed = True)
134       or not (closed = False)
135
136    Outputs:
137       indices: array of same length as points with indices of points falling
138       inside the polygon listed from the beginning and indices of points
139       falling outside listed from the end.
140
141       count: count of points falling inside the polygon
142
143       The indices of points inside are obtained as indices[:count]
144       The indices of points outside are obtained as indices[count:]
145
146
147    Examples:
148       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
149
150       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
151       will return the indices [0, 2, 1] and count == 2 as only the first
152       and the last point are inside the unit square
153
154    Remarks:
155       The vertices may be listed clockwise or counterclockwise and
156       the first point may optionally be repeated.
157       Polygons do not need to be convex.
158       Polygons can have holes in them and points inside a hole is
159       regarded as being outside the polygon.
160
161    Algorithm is based on work by Darel Finley,
162    http://www.alienryderflex.com/polygon/
163
164    """
165
166    from Numeric import array, Float, reshape, Int, zeros
167
168
169    #Input checks
170    try:
171        points = ensure_numeric(points, Float)
172    except:
173        msg = 'Points could not be converted to Numeric array'
174        raise msg
175
176    try:
177        polygon = ensure_numeric(polygon, Float)
178    except:
179        msg = 'Polygon could not be converted to Numeric array'
180        raise msg
181
182    assert len(polygon.shape) == 2,\
183       'Polygon array must be a 2d array of vertices'
184
185    assert polygon.shape[1] == 2,\
186       'Polygon array must have two columns'
187
188    assert len(points.shape) == 2,\
189       'Points array must be a 2d array'
190
191    assert points.shape[1] == 2,\
192       'Points array must have two columns'
193
194    N = polygon.shape[0] #Number of vertices in polygon
195    M = points.shape[0]  #Number of points
196
197    px = polygon[:,0]
198    py = polygon[:,1]
199
200    #Used for an optimisation when points are far away from polygon
201    minpx = min(px); maxpx = max(px)
202    minpy = min(py); maxpy = max(py)
203
204
205    #Begin main loop
206    indices = zeros(M, Int)
207
208    inside_index = 0    #Keep track of points inside
209    outside_index = M-1 #Keep track of points outside (starting from end)
210
211    for k in range(M):
212
213        if verbose:
214            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
215
216        #for each point
217        x = points[k, 0]
218        y = points[k, 1]
219
220        inside = False
221
222        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
223            #Check polygon
224            for i in range(N):
225                j = (i+1)%N
226
227                #Check for case where point is contained in line segment
228                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
229                    if closed:
230                        inside = True
231                    else:
232                        inside = False
233                    break
234                else:
235                    #Check if truly inside polygon
236                    if py[i] < y and py[j] >= y or\
237                       py[j] < y and py[i] >= y:
238                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
239                            inside = not inside
240
241        if inside:
242            indices[inside_index] = k
243            inside_index += 1
244        else:
245            indices[outside_index] = k
246            outside_index -= 1
247
248    return indices, inside_index
249
250
251def separate_points_by_polygon_c(points, polygon,
252                                 closed = True, verbose = False):
253    """Determine whether points are inside or outside a polygon
254
255    C-wrapper
256    """
257
258    from Numeric import array, Float, reshape, zeros, Int
259
260
261    if verbose: print 'Checking input to separate_points_by_polygon'
262    #Input checks
263    try:
264        points = ensure_numeric(points, Float)
265    except:
266        msg = 'Points could not be converted to Numeric array'
267        raise msg
268
269    #if verbose: print 'Checking input to separate_points_by_polygon 2'
270    try:
271        polygon = ensure_numeric(polygon, Float)
272    except:
273        msg = 'Polygon could not be converted to Numeric array'
274        raise msg
275
276    if verbose: print 'check'
277
278    assert len(polygon.shape) == 2,\
279       'Polygon array must be a 2d array of vertices'
280
281    assert polygon.shape[1] == 2,\
282       'Polygon array must have two columns'
283
284    assert len(points.shape) == 2,\
285       'Points array must be a 2d array'
286
287    assert points.shape[1] == 2,\
288       'Points array must have two columns'
289
290    N = polygon.shape[0] #Number of vertices in polygon
291    M = points.shape[0]  #Number of points
292
293    from polygon_ext import separate_points_by_polygon
294
295    if verbose: print 'Allocating array for indices'
296
297    indices = zeros( M, Int )
298
299    if verbose: print 'Calling C-version of inside poly'
300    count = separate_points_by_polygon(points, polygon, indices,
301                                       int(closed), int(verbose))
302
303    return indices, count
304
305
306
307class Polygon_function:
308    """Create callable object f: x,y -> z, where a,y,z are vectors and
309    where f will return different values depending on whether x,y belongs
310    to specified polygons.
311
312    To instantiate:
313
314       Polygon_function(polygons)
315
316    where polygons is a list of tuples of the form
317
318      [ (P0, v0), (P1, v1), ...]
319
320      with Pi being lists of vertices defining polygons and vi either
321      constants or functions of x,y to be applied to points with the polygon.
322
323    The function takes an optional argument, default which is the value
324    (or function) to used for points not belonging to any polygon.
325    For example:
326
327       Polygon_function(polygons, default = 0.03)
328
329    If omitted the default value will be 0.0
330
331    Note: If two polygons overlap, the one last in the list takes precedence
332
333    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???
334
335    """
336
337    def __init__(self, regions, default = 0.0):
338
339        try:
340            len(regions)
341        except:
342            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
343            raise msg
344
345
346        T = regions[0]
347        try:
348            a = len(T)
349        except:
350            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
351            raise msg
352
353        assert a == 2, 'Must have two component each: %s' %T
354
355        self.regions = regions
356        self.default = default
357
358
359    def __call__(self, x, y):
360        from Numeric import ones, Float, concatenate, array, reshape, choose
361
362        x = array(x).astype(Float)
363        y = array(y).astype(Float)
364
365        N = len(x)
366        assert len(y) == N
367
368        points = concatenate( (reshape(x, (N, 1)),
369                               reshape(y, (N, 1))), axis=1 )
370
371        if callable(self.default):
372            z = self.default(x,y)
373        else:
374            z = ones(N, Float) * self.default
375
376        for polygon, value in self.regions:
377            indices = inside_polygon(points, polygon)
378
379            #FIXME: This needs to be vectorised
380            if callable(value):
381                for i in indices:
382                    xx = array([x[i]])
383                    yy = array([y[i]])
384                    z[i] = value(xx, yy)[0]
385            else:
386                for i in indices:
387                    z[i] = value
388
389        return z
390
392    """Read points assumed to form a polygon
393       There must be exactly two numbers in each line
394    """
395
396    #Get polygon
397    fid = open(filename)
399    fid.close()
400    polygon = []
401    for line in lines:
402        fields = line.split(',')
403        polygon.append( [float(fields[0]), float(fields[1])] )
404
405    return polygon
406
407def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
408    """Populate given polygon with uniformly distributed points.
409
410    Input:
411       polygon - list of vertices of polygon
412       number_of_points - (optional) number of points
413       seed - seed for random number generator (default=None)
414       exclude - list of polygons (inside main polygon) from where points should be excluded
415
416    Output:
417       points - list of points inside polygon
418
419    Examples:
420       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
421       will return five randomly selected points inside the unit square
422    """
423
424    from random import uniform, seed as seed_function
425
426    seed_function(seed)
427
428    points = []
429
430    #Find outer extent of polygon
431    max_x = min_x = polygon[0][0]
432    max_y = min_y = polygon[0][1]
433    for point in polygon[1:]:
434        x = point[0]
435        if x > max_x: max_x = x
436        if x < min_x: min_x = x
437        y = point[1]
438        if y > max_y: max_y = y
439        if y < min_y: min_y = y
440
441
442    while len(points) < number_of_points:
443        x = uniform(min_x, max_x)
444        y = uniform(min_y, max_y)
445
446        append = False
447        if inside_polygon( [x,y], polygon ):
448
449            append = True
450
451            #Check exclusions
452            if exclude is not None:
453                for ex_poly in exclude:
454                    if inside_polygon( [x,y], ex_poly ):
455                        append = False
456
457
458        if append is True:
459            points.append([x,y])
460
461    return points
462
463def point_in_polygon(polygon, delta=1e-8):
464    """Return a point inside a given polygon which will be close to the
465    polygon edge.
466
467    Input:
468       polygon - list of vertices of polygon
469       delta - the square root of 2 * delta is the maximum distance from the
470       polygon points and the returned point.
471    Output:
472       points - a point inside polygon
473
474       searches in all diagonals and up and down (not left and right)
475    """
476    import exceptions
477    class Found(exceptions.Exception): pass
478
479    point_in = False
480    while not point_in:
481        try:
482            for poly_point in polygon: #[1:]:
483                for x_mult in range (-1,2):
484                    for y_mult in range (-1,2):
485                        x = poly_point[0]
486                        y = poly_point[1]
487                        if x == 0:
488                            x_delta = x_mult*delta
489                        else:
490                            x_delta = x+x_mult*x*delta
491
492                        if y == 0:
493                            y_delta = y_mult*delta
494                        else:
495                            y_delta = y+y_mult*y*delta
496
497                        point = [x_delta, y_delta]
498                        #print "point",point
499                        if inside_polygon( point, polygon, closed = False ):
500                            raise Found
501        except Found:
502            point_in = True
503        else:
504            delta = delta*0.1
505    return point
506
507
508
509##############################################
510#Initialise module
511
512import compile
513if compile.can_use_C_extension('polygon_ext.c'):
514    from polygon_ext import point_on_line
515    separate_points_by_polygon = separate_points_by_polygon_c
516
517
518if __name__ == "__main__":
519    pass
Note: See TracBrowser for help on using the repository browser.