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

Last change on this file since 2544 was 2531, checked in by ole, 19 years ago

Investigating feasibility of scipy as an alternative to Numeric

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