# source:anuga_core/source/anuga/utilities/polygon.py@3633

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

Fixed diagnostic counting error and improved messages.

File size: 20.7 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 anuga.utilities.numerical_tools import ensure_numeric
16from anuga.geospatial_data.geospatial_data import ensure_absolute
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, -a])
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 is_inside_polygon(point, polygon, closed=True, verbose=False):
46    """Determine if one point is inside a polygon
47
48    See inside_polygon for more details
49    """
50
51    indices = inside_polygon(point, polygon, closed, verbose)
52
53    if indices.shape == 1:
54        return True
55    elif indices.shape == 0:
56        return False
57    else:
58        msg = 'is_inside_polygon must be invoked with one point only'
59        raise msg
60
61
62def inside_polygon(points, polygon, closed=True, verbose=False):
63    """Determine points inside a polygon
64
65       Functions inside_polygon and outside_polygon have been defined in
66       terms af separate_by_polygon which will put all inside indices in
67       the first part of the indices array and outside indices in the last
68
69       See separate_points_by_polygon for documentation
70
71       points and polygon can be a geospatial instance,
72       a list or a numeric array
73    """
74
75    if verbose: print 'Checking input to inside_polygon'
76
77    try:
78        points = ensure_absolute(points)
79    except NameError, e:
80        raise NameError, e
81    except:
82        msg = 'Points could not be converted to Numeric array'
83        raise msg
84
85    try:
86        polygon = ensure_absolute(polygon)
87    except NameError, e:
88        raise NameError, e
89    except:
90        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
91        raise msg
92
93    if len(points.shape) == 1:
94        # Only one point was passed in. Convert to array of points
95        points = reshape(points, (1,2))
96
97    indices, count = separate_points_by_polygon(points, polygon,
98                                                closed=closed,
99                                                verbose=verbose)
100
101    # Return indices of points inside polygon
102    return indices[:count]
103
104
105
106def is_outside_polygon(point, polygon, closed=True, verbose=False,
107                       points_geo_ref=None, polygon_geo_ref=None):
108    """Determine if one point is outside a polygon
109
110    See outside_polygon for more details
111    """
112
113    indices = outside_polygon(point, polygon, closed, verbose)
114                              #points_geo_ref, polygon_geo_ref)
115
116    if indices.shape == 1:
117        return True
118    elif indices.shape == 0:
119        return False
120    else:
121        msg = 'is_outside_polygon must be invoked with one point only'
122        raise msg
123
124
125def outside_polygon(points, polygon, closed = True, verbose = False):
126    """Determine points outside a polygon
127
128       Functions inside_polygon and outside_polygon have been defined in
129       terms af separate_by_polygon which will put all inside indices in
130       the first part of the indices array and outside indices in the last
131
132       See separate_points_by_polygon for documentation
133    """
134
135    if verbose: print 'Checking input to outside_polygon'
136    try:
137        points = ensure_numeric(points, Float)
138    except NameError, e:
139        raise NameError, e
140    except:
141        msg = 'Points could not be converted to Numeric array'
142        raise msg
143
144    try:
145        polygon = ensure_numeric(polygon, Float)
146    except NameError, e:
147        raise NameError, e
148    except:
149        msg = 'Polygon could not be converted to Numeric array'
150        raise msg
151
152
153    if len(points.shape) == 1:
154        # Only one point was passed in. Convert to array of points
155        points = reshape(points, (1,2))
156
157    indices, count = separate_points_by_polygon(points, polygon,
158                                                closed=closed,
159                                                verbose=verbose)
160
161    # Return indices of points outside polygon
162    if count == len(indices):
163        # No points are outside
164        return array([])
165    else:
166        return indices[count:][::-1]  #return reversed
167
168
169def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
170    """Determine points inside and outside a polygon
171
172       See separate_points_by_polygon for documentation
173
174       Returns an array of points inside and an array of points outside the polygon
175    """
176
177    if verbose: print 'Checking input to outside_polygon'
178    try:
179        points = ensure_numeric(points, Float)
180    except NameError, e:
181        raise NameError, e
182    except:
183        msg = 'Points could not be converted to Numeric array'
184        raise msg
185
186    try:
187        polygon = ensure_numeric(polygon, Float)
188    except NameError, e:
189        raise NameError, e
190    except:
191        msg = 'Polygon could not be converted to Numeric array'
192        raise msg
193
194    if len(points.shape) == 1:
195        # Only one point was passed in. Convert to array of points
196        points = reshape(points, (1,2))
197
198
199    indices, count = separate_points_by_polygon(points, polygon,
200                                                closed=closed,
201                                                verbose=verbose)
202
203    # Returns indices of points inside and indices of points outside
204    # the polygon
205
206    if count == len(indices):
207        # No points are outside
208        return indices[:count],[]
209    else:
210        return  indices[:count], indices[count:][::-1]  #return reversed
211
212
213def separate_points_by_polygon(points, polygon,
214                               closed = True, verbose = False):
215    """Determine whether points are inside or outside a polygon
216
217    Input:
218       points - Tuple of (x, y) coordinates, or list of tuples
219       polygon - list of vertices of polygon
220       closed - (optional) determine whether points on boundary should be
221       regarded as belonging to the polygon (closed = True)
222       or not (closed = False)
223
224    Outputs:
225       indices: array of same length as points with indices of points falling
226       inside the polygon listed from the beginning and indices of points
227       falling outside listed from the end.
228
229       count: count of points falling inside the polygon
230
231       The indices of points inside are obtained as indices[:count]
232       The indices of points outside are obtained as indices[count:]
233
234
235    Examples:
236       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
237
238       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
239       will return the indices [0, 2, 1] and count == 2 as only the first
240       and the last point are inside the unit square
241
242    Remarks:
243       The vertices may be listed clockwise or counterclockwise and
244       the first point may optionally be repeated.
245       Polygons do not need to be convex.
246       Polygons can have holes in them and points inside a hole is
247       regarded as being outside the polygon.
248
249    Algorithm is based on work by Darel Finley,
250    http://www.alienryderflex.com/polygon/
251
252    """
253
254    #Input checks
255
256
257    try:
258        points = ensure_numeric(points, Float)
259    except NameError, e:
260        raise NameError, e
261    except:
262        msg = 'Points could not be converted to Numeric array'
263        raise msg
264
265    try:
266        polygon = ensure_numeric(polygon, Float)
267    except NameError, e:
268        raise NameError, e
269    except:
270        msg = 'Polygon could not be converted to Numeric array'
271        raise msg
272
273    assert len(polygon.shape) == 2,\
274       'Polygon array must be a 2d array of vertices'
275
276    assert polygon.shape == 2,\
277       'Polygon array must have two columns'
278
279    assert len(points.shape) == 2,\
280       'Points array must be a 2d array'
281
282    assert points.shape == 2,\
283       'Points array must have two columns'
284
285    N = polygon.shape #Number of vertices in polygon
286    M = points.shape  #Number of points
287
288    px = polygon[:,0]
289    py = polygon[:,1]
290
291    #Used for an optimisation when points are far away from polygon
292    minpx = min(px); maxpx = max(px)
293    minpy = min(py); maxpy = max(py)
294
295
296    #Begin main loop
297    indices = zeros(M, Int)
298
299    inside_index = 0    #Keep track of points inside
300    outside_index = M-1 #Keep track of points outside (starting from end)
301
302    if verbose: print 'Separating %d points' %M
303    for k in range(M):
304
305        if verbose:
306            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
307
308        #for each point
309        x = points[k, 0]
310        y = points[k, 1]
311
312        inside = False
313
314        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
315            #Check polygon
316            for i in range(N):
317                j = (i+1)%N
318
319                #Check for case where point is contained in line segment
320                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
321                    if closed:
322                        inside = True
323                    else:
324                        inside = False
325                    break
326                else:
327                    #Check if truly inside polygon
328                    if py[i] < y and py[j] >= y or\
329                       py[j] < y and py[i] >= y:
330                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
331                            inside = not inside
332
333        if inside:
334            indices[inside_index] = k
335            inside_index += 1
336        else:
337            indices[outside_index] = k
338            outside_index -= 1
339
340    return indices, inside_index
341
342
343def separate_points_by_polygon_c(points, polygon,
344                                 closed = True, verbose = False):
345    """Determine whether points are inside or outside a polygon
346
347    C-wrapper
348    """
349
350
351    if verbose: print 'Checking input to separate_points_by_polygon'
352
353
354    #Input checks
355
356    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
357    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
358
359
360    try:
361        points = ensure_numeric(points, Float)
362    except NameError, e:
363        raise NameError, e
364    except:
365        msg = 'Points could not be converted to Numeric array'
366        raise msg
367
368    #if verbose: print 'Checking input to separate_points_by_polygon 2'
369    try:
370        polygon = ensure_numeric(polygon, Float)
371    except NameError, e:
372        raise NameError, e
373    except:
374        msg = 'Polygon could not be converted to Numeric array'
375        raise msg
376
377    #if verbose: print 'check'
378
379    assert len(polygon.shape) == 2,\
380       'Polygon array must be a 2d array of vertices'
381
382    assert polygon.shape == 2,\
383       'Polygon array must have two columns'
384
385    assert len(points.shape) == 2,\
386       'Points array must be a 2d array'
387
388    assert points.shape == 2,\
389       'Points array must have two columns'
390
391    N = polygon.shape #Number of vertices in polygon
392    M = points.shape  #Number of points
393
394    from polygon_ext import separate_points_by_polygon
395
396    if verbose: print 'Allocating array for indices'
397
398    indices = zeros( M, Int )
399
400    #if verbose: print 'Calling C-version of inside poly'
401
402    count = separate_points_by_polygon(points, polygon, indices,
403                                       int(closed), int(verbose))
404
405    if verbose: print 'Found %d points (out of %d) inside polygon'\
406       %(count, M)
407    return indices, count
408
409
410def polygon_area(polygon):
411    """ Determin area of arbitrary polygon
412    Reference
413    http://mathworld.wolfram.com/PolygonArea.html
414    """
415
416    n = len(polygon)
417    poly_area = 0.0
418
419    for i in range(n):
420        pti = polygon[i]
421        if i == n-1:
422            pt1 = polygon
423        else:
424            pt1 = polygon[i+1]
425        xi = pti
426        yi1 = pt1
427        xi1 = pt1
428        yi = pti
429        poly_area += xi*yi1 - xi1*yi
430
431    return abs(poly_area/2)
432
433def plot_polygons(polygons, figname, verbose=False):
434
435    """ Take list of polygons and plot.
436
437    Inputs:
438
439    polygons         - list of polygons
440
441    figname          - name to save figure to
442
443    Outputs:
444
445    - list of min and max of x and y coordinates
446    - plot of polygons
447    """
448
449    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close
450
451    assert type(polygons) == list,\
452               'input must be a list of polygons'
453
454    ion()
455    hold(True)
456
457    minx = 1e10
458    maxx = 0.0
459    miny = 1e10
460    maxy = 0.0
461
462    for polygon in polygons:
463        x, y = poly_xy(polygon)
464        if min(x) < minx: minx = min(x)
465        if max(x) > maxx: maxx = max(x)
466        if min(y) < miny: miny = min(y)
467        if max(y) > maxy: maxy = max(y)
468        plot(x,y,'r-')
469        xlabel('x')
470        ylabel('y')
471
472    savefig(figname)
473
474    close('all')
475
476    vec = [minx,maxx,miny,maxy]
477
478    return vec
479
480def poly_xy(polygon, verbose=False):
481    """ this is used within plot_polygons so need to duplicate
482        the first point so can have closed polygon in plot
483    """
484
485    if verbose: print 'Checking input to poly_xy'
486
487    try:
488        polygon = ensure_numeric(polygon, Float)
489    except NameError, e:
490        raise NameError, e
491    except:
492        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
493        raise msg
494
495    x = polygon[:,0]
496    y = polygon[:,1]
497    x = concatenate((x, [polygon[0,0]]), axis = 0)
498    y = concatenate((y, [polygon[0,1]]), axis = 0)
499
500    return x, y
501
502#    x = []
503#    y = []
504#    n = len(poly)
505#    firstpt = poly
506#    for i in range(n):
507#        thispt = poly[i]
508#        x.append(thispt)
509#        y.append(thispt)
510
511#    x.append(firstpt)
512#    y.append(firstpt)
513
514#    return x, y
515
516class Polygon_function:
517    """Create callable object f: x,y -> z, where a,y,z are vectors and
518    where f will return different values depending on whether x,y belongs
519    to specified polygons.
520
521    To instantiate:
522
523       Polygon_function(polygons)
524
525    where polygons is a list of tuples of the form
526
527      [ (P0, v0), (P1, v1), ...]
528
529      with Pi being lists of vertices defining polygons and vi either
530      constants or functions of x,y to be applied to points with the polygon.
531
532    The function takes an optional argument, default which is the value
533    (or function) to used for points not belonging to any polygon.
534    For example:
535
536       Polygon_function(polygons, default = 0.03)
537
538    If omitted the default value will be 0.0
539
540    Note: If two polygons overlap, the one last in the list takes precedence
541
542    Coordinates specified in the call are assumed to be relative to the origin (georeference)
543    e.g. used by domain. By specifying the optional argument georeference, all points are made relative.
544
545    FIXME: This should really work with geo_spatial point sets.
546    """
547
548    def __init__(self, regions, default = 0.0, geo_reference = None):
549
550        try:
551            len(regions)
552        except:
553            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
554            raise msg
555
556
557        T = regions
558        try:
559            a = len(T)
560        except:
561            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
562            raise msg
563
564        assert a == 2, 'Must have two component each: %s' %T
565
566
567        if geo_reference is None:
568            from anuga.coordinate_transforms.geo_reference import Geo_reference
569            geo_reference = Geo_reference()
570
571
572        self.default = default
573
574        #Make points in polygons relative to geo_reference
575        self.regions = []
576        for polygon, value in regions:
577            P = geo_reference.change_points_geo_ref(polygon)
578            self.regions.append( (P, value) )
579
580
581
582
583    def __call__(self, x, y):
584        x = array(x).astype(Float)
585        y = array(y).astype(Float)
586
587        N = len(x)
588        assert len(y) == N
589
590        points = concatenate( (reshape(x, (N, 1)),
591                               reshape(y, (N, 1))), axis=1 )
592
593        if callable(self.default):
594            z = self.default(x,y)
595        else:
596            z = ones(N, Float) * self.default
597
598        for polygon, value in self.regions:
599            indices = inside_polygon(points, polygon)
600
601            #FIXME: This needs to be vectorised
602            if callable(value):
603                for i in indices:
604                    xx = array([x[i]])
605                    yy = array([y[i]])
606                    z[i] = value(xx, yy)
607            else:
608                for i in indices:
609                    z[i] = value
610
611        return z
612
614    """Read points assumed to form a polygon.
615       There must be exactly two numbers in each line separated by a comma.
617    """
618
619    #Get polygon
620    fid = open(filename)
622    fid.close()
623    polygon = []
624    for line in lines:
625        fields = line.split(split)
626        polygon.append( [float(fields), float(fields)] )
627
628    return polygon
629
630def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
631    """Populate given polygon with uniformly distributed points.
632
633    Input:
634       polygon - list of vertices of polygon
635       number_of_points - (optional) number of points
636       seed - seed for random number generator (default=None)
637       exclude - list of polygons (inside main polygon) from where points should be excluded
638
639    Output:
640       points - list of points inside polygon
641
642    Examples:
643       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
644       will return five randomly selected points inside the unit square
645    """
646
647    from random import uniform, seed as seed_function
648
649    seed_function(seed)
650
651    points = []
652
653    #Find outer extent of polygon
654    max_x = min_x = polygon
655    max_y = min_y = polygon
656    for point in polygon[1:]:
657        x = point
658        if x > max_x: max_x = x
659        if x < min_x: min_x = x
660        y = point
661        if y > max_y: max_y = y
662        if y < min_y: min_y = y
663
664
665    while len(points) < number_of_points:
666        x = uniform(min_x, max_x)
667        y = uniform(min_y, max_y)
668
669        append = False
670        if is_inside_polygon([x,y], polygon):
671
672            append = True
673
674            #Check exclusions
675            if exclude is not None:
676                for ex_poly in exclude:
677                    if is_inside_polygon([x,y], ex_poly):
678                        append = False
679
680
681        if append is True:
682            points.append([x,y])
683
684    return points
685
686def point_in_polygon(polygon, delta=1e-8):
687    """Return a point inside a given polygon which will be close to the
688    polygon edge.
689
690    Input:
691       polygon - list of vertices of polygon
692       delta - the square root of 2 * delta is the maximum distance from the
693       polygon points and the returned point.
694    Output:
695       points - a point inside polygon
696
697       searches in all diagonals and up and down (not left and right)
698    """
699    import exceptions
700    class Found(exceptions.Exception): pass
701
702    point_in = False
703    while not point_in:
704        try:
705            for poly_point in polygon: #[1:]:
706                for x_mult in range (-1,2):
707                    for y_mult in range (-1,2):
708                        x = poly_point
709                        y = poly_point
710                        if x == 0:
711                            x_delta = x_mult*delta
712                        else:
713                            x_delta = x+x_mult*x*delta
714
715                        if y == 0:
716                            y_delta = y_mult*delta
717                        else:
718                            y_delta = y+y_mult*y*delta
719
720                        point = [x_delta, y_delta]
721                        #print "point",point
722                        if is_inside_polygon(point, polygon, closed=False):
723                            raise Found
724        except Found:
725            point_in = True
726        else:
727            delta = delta*0.1
728    return point
729
730
731
732##############################################
733#Initialise module
734
735import compile
736if compile.can_use_C_extension('polygon_ext.c'):
737    from polygon_ext import point_on_line
738    separate_points_by_polygon = separate_points_by_polygon_c
739
740
741if __name__ == "__main__":
742    pass
Note: See TracBrowser for help on using the repository browser.