source: anuga_core/source/anuga/utilities/polygon.py @ 4669

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

Typo

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