source: inundation-numpy-branch/utilities/polygon.py @ 3381

Last change on this file since 3381 was 2548, checked in by ole, 19 years ago

A few more fixes towards numpy

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