source: inundation/utilities/polygon.py @ 2375

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

Added georeferencing to polygon_function and added test

File size: 14.8 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 %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[1] == 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[1] == 2,\
194       'Points array must have two columns'
195
196    N = polygon.shape[0] #Number of vertices in polygon
197    M = points.shape[0]  #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[1] == 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[1] == 2,\
290       'Points array must have two columns'
291
292    N = polygon.shape[0] #Number of vertices in polygon
293    M = points.shape[0]  #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    Coordinates specified in the call are assumed to be relative to the origin (georeference)
336    e.g. used by domain. By specifying the optional argument georeference, all points are made relative.
337
338    FIXME: This should really work with geo_spatial point sets.
339    """
340
341    def __init__(self, regions, default = 0.0, geo_reference = None):
342
343        try:
344            len(regions)
345        except:
346            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
347            raise msg
348
349
350        T = regions[0]
351        try:
352            a = len(T)
353        except:
354            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
355            raise msg
356
357        assert a == 2, 'Must have two component each: %s' %T
358
359
360        if geo_reference is None:
361            from coordinate_transforms.geo_reference import Geo_reference
362            geo_reference = Geo_reference()
363
364
365        self.default = default
366
367        #Make points in polygons relative to geo_reference
368        self.regions = []           
369        for polygon, value in regions:
370            P = geo_reference.change_points_geo_ref(polygon)           
371            self.regions.append( (P, value) )
372
373
374
375
376    def __call__(self, x, y):
377        from Numeric import ones, Float, concatenate, array, reshape, choose
378
379        x = array(x).astype(Float)
380        y = array(y).astype(Float)
381
382        N = len(x)
383        assert len(y) == N
384
385        points = concatenate( (reshape(x, (N, 1)),
386                               reshape(y, (N, 1))), axis=1 )
387
388        if callable(self.default):
389            z = self.default(x,y)
390        else:
391            z = ones(N, Float) * self.default
392
393        for polygon, value in self.regions:
394            indices = inside_polygon(points, polygon)
395
396            #FIXME: This needs to be vectorised
397            if callable(value):
398                for i in indices:
399                    xx = array([x[i]])
400                    yy = array([y[i]])
401                    z[i] = value(xx, yy)[0]
402            else:
403                for i in indices:
404                    z[i] = value
405
406        return z
407
408def read_polygon(filename):
409    """Read points assumed to form a polygon.
410       There must be exactly two numbers in each line separated by a comma.
411       No header.
412    """
413
414    #Get polygon
415    fid = open(filename)
416    lines = fid.readlines()
417    fid.close()
418    polygon = []
419    for line in lines:
420        fields = line.split(',')
421        polygon.append( [float(fields[0]), float(fields[1])] )
422
423    return polygon
424
425def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
426    """Populate given polygon with uniformly distributed points.
427
428    Input:
429       polygon - list of vertices of polygon
430       number_of_points - (optional) number of points
431       seed - seed for random number generator (default=None)
432       exclude - list of polygons (inside main polygon) from where points should be excluded
433
434    Output:
435       points - list of points inside polygon
436
437    Examples:
438       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
439       will return five randomly selected points inside the unit square
440    """
441
442    from random import uniform, seed as seed_function
443
444    seed_function(seed)
445
446    points = []
447
448    #Find outer extent of polygon
449    max_x = min_x = polygon[0][0]
450    max_y = min_y = polygon[0][1]
451    for point in polygon[1:]:
452        x = point[0]
453        if x > max_x: max_x = x
454        if x < min_x: min_x = x
455        y = point[1]
456        if y > max_y: max_y = y
457        if y < min_y: min_y = y
458
459
460    while len(points) < number_of_points:
461        x = uniform(min_x, max_x)
462        y = uniform(min_y, max_y)
463
464        append = False
465        if inside_polygon( [x,y], polygon ):
466
467            append = True
468           
469            #Check exclusions
470            if exclude is not None:
471                for ex_poly in exclude:
472                    if inside_polygon( [x,y], ex_poly ):
473                        append = False
474                       
475                   
476        if append is True:               
477            points.append([x,y])               
478
479    return points
480
481def point_in_polygon(polygon, delta=1e-8):
482    """Return a point inside a given polygon which will be close to the
483    polygon edge.
484
485    Input:
486       polygon - list of vertices of polygon
487       delta - the square root of 2 * delta is the maximum distance from the
488       polygon points and the returned point.
489    Output:
490       points - a point inside polygon
491
492       searches in all diagonals and up and down (not left and right)
493    """
494    import exceptions
495    class Found(exceptions.Exception): pass
496
497    point_in = False
498    while not point_in:
499        try:
500            for poly_point in polygon: #[1:]:
501                for x_mult in range (-1,2):
502                    for y_mult in range (-1,2):
503                        x = poly_point[0]
504                        y = poly_point[1]
505                        if x == 0:
506                            x_delta = x_mult*delta
507                        else:
508                            x_delta = x+x_mult*x*delta
509                           
510                        if y == 0:
511                            y_delta = y_mult*delta
512                        else:
513                            y_delta = y+y_mult*y*delta
514                           
515                        point = [x_delta, y_delta]
516                        #print "point",point
517                        if inside_polygon( point, polygon, closed = False ):
518                            raise Found
519        except Found: 
520            point_in = True
521        else:
522            delta = delta*0.1 
523    return point
524
525
526
527##############################################
528#Initialise module
529
530import compile
531if compile.can_use_C_extension('polygon_ext.c'):
532    from polygon_ext import point_on_line
533    separate_points_by_polygon = separate_points_by_polygon_c
534
535
536if __name__ == "__main__":
537    pass
Note: See TracBrowser for help on using the repository browser.