# source:inundation/utilities/polygon.py@2351

Last change on this file since 2351 was 2351, checked in by ole, 17 years ago

Minor bug fix

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