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

Last change on this file since 4768 was 4768, checked in by ole, 16 years ago

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