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

Last change on this file since 5225 was 5225, checked in by ole, 15 years ago

Work done during Water Down Under 2008.
Line Mesh intersections towards ticket:175 (flow through a cross section).

File size: 23.9 KB
Line
1#!/usr/bin/env python
2"""Polygon manipulations
3
4"""
5
6
7#try:
8#    from scipy import Float, Int, zeros, ones, array, concatenate, reshape, dot
9#except:
10#    #print 'Could not find scipy - using Numeric'
11
12from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose
13
14
15from math import sqrt
16from anuga.utilities.numerical_tools import ensure_numeric
17from anuga.geospatial_data.geospatial_data import ensure_absolute
18
19def point_on_line(point, line):
20    """Determine whether a point is on a line segment
21
22    Input:
23        point is given by [x, y]
24        line is given by [x0, y0], [x1, y1]] or
25        the equivalent 2x2 Numeric array with each row corresponding to a point.
26
27    Output:
28
29    """
30
31    point = ensure_numeric(point)
32    line = ensure_numeric(line)
33
34
35    res = _point_on_line(point[0], point[1],
36                         line[0,0], line[0,1],
37                         line[1,0], line[1,1])
38
39    return bool(res)
40
41
42
43
44
45def intersection(line0, line1):
46    """Returns intersecting point between two line segments or None (if parallel or no intersection is found)
47
48    Inputs:
49        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
50                      A line can also be a 2x2 numeric array with each row corresponding to a point.
51
52
53    Output:
54        Point [x,y] or None.
55
56    If line extensions intersect outside their limits, None is returned as well.
57
58    """
59
60    # FIXME (Ole): Write this in C
61
62    line0 = ensure_numeric(line0, Float)
63    line1 = ensure_numeric(line1, Float)
64
65    x0 = line0[0,0]; y0 = line0[0,1]
66    x1 = line0[1,0]; y1 = line0[1,1]
67
68    x2 = line1[0,0]; y2 = line1[0,1]
69    x3 = line1[1,0]; y3 = line1[1,1]
70
71    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
72
73    #print 'denom', denom, line0, line1
74    if allclose(denom, 0.0):
75        # Lines are parallel
76        return None
77    else:
78        u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
79        u0 = u0/denom
80
81        u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
82        u1 = u1/denom
83
84        x = x0 + u0*(x1-x0)
85        y = y0 + u0*(y1-y0)
86
87        assert allclose(x, x2 + u1*(x3-x2))
88        assert allclose(y, y2 + u1*(y3-y2))
89
90        # Check if point found lies within given line segments
91        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
92            # We have intersection
93
94            # Need tolerances if going ahead with this check
95            #msg = 'Algorithm error. Intersection was detected but point (%f, %f) does not lie ' %(x,y)
96            #msg += 'on line0: %s' %(line0)
97            #assert point_on_line([x, y], line0), msg
98            #msg = 'Algorithm error. Intersection was detected but point (%f, %f) does not lie ' %(x,y)
99            #msg += 'on line1: %s' %(line1)
100            #assert point_on_line([x, y], line1), msg
101
102            return [x, y]
103        else:
104            return None
105
106
107def is_inside_polygon(point, polygon, closed=True, verbose=False):
108    """Determine if one point is inside a polygon
109
110    See inside_polygon for more details
111    """
112
113    indices = inside_polygon(point, polygon, closed, verbose)
114
115    if indices.shape[0] == 1:
116        return True
117    elif indices.shape[0] == 0:
118        return False
119    else:
120        msg = 'is_inside_polygon must be invoked with one point only'
121        raise msg
122
123
124def inside_polygon(points, polygon, closed=True, verbose=False):
125    """Determine points inside a polygon
126
127       Functions inside_polygon and outside_polygon have been defined in
128       terms af separate_by_polygon which will put all inside indices in
129       the first part of the indices array and outside indices in the last
130
131       See separate_points_by_polygon for documentation
132
133       points and polygon can be a geospatial instance,
134       a list or a numeric array
135    """
136
137    #if verbose: print 'Checking input to inside_polygon'
138
139    try:
140        points = ensure_absolute(points)
141    except NameError, e:
142        raise NameError, e
143    except:
144        # If this fails it is going to be because the points can't be
145        # converted to a numeric array.
146        msg = 'Points could not be converted to Numeric array'
147        raise msg
148
149    try:
150        polygon = ensure_absolute(polygon)
151    except NameError, e:
152        raise NameError, e
153    except:
154        # If this fails it is going to be because the points can't be
155        # converted to a numeric array.
156        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
157        raise msg
158
159    if len(points.shape) == 1:
160        # Only one point was passed in. Convert to array of points
161        points = reshape(points, (1,2))
162
163    indices, count = separate_points_by_polygon(points, polygon,
164                                                closed=closed,
165                                                verbose=verbose)
166
167    # Return indices of points inside polygon
168    return indices[:count]
169
170
171
172def is_outside_polygon(point, polygon, closed=True, verbose=False,
173                       points_geo_ref=None, polygon_geo_ref=None):
174    """Determine if one point is outside a polygon
175
176    See outside_polygon for more details
177    """
178
179    indices = outside_polygon(point, polygon, closed, verbose)
180                              #points_geo_ref, polygon_geo_ref)
181
182    if indices.shape[0] == 1:
183        return True
184    elif indices.shape[0] == 0:
185        return False
186    else:
187        msg = 'is_outside_polygon must be invoked with one point only'
188        raise msg
189
190
191def outside_polygon(points, polygon, closed = True, verbose = False):
192    """Determine points outside a polygon
193
194       Functions inside_polygon and outside_polygon have been defined in
195       terms af separate_by_polygon which will put all inside indices in
196       the first part of the indices array and outside indices in the last
197
198       See separate_points_by_polygon for documentation
199    """
200
201    #if verbose: print 'Checking input to outside_polygon'
202    try:
203        points = ensure_numeric(points, Float)
204    except NameError, e:
205        raise NameError, e
206    except:
207        msg = 'Points could not be converted to Numeric array'
208        raise msg
209
210    try:
211        polygon = ensure_numeric(polygon, Float)
212    except NameError, e:
213        raise NameError, e
214    except:
215        msg = 'Polygon could not be converted to Numeric array'
216        raise msg
217
218
219    if len(points.shape) == 1:
220        # Only one point was passed in. Convert to array of points
221        points = reshape(points, (1,2))
222
223    indices, count = separate_points_by_polygon(points, polygon,
224                                                closed=closed,
225                                                verbose=verbose)
226
227    # Return indices of points outside polygon
228    if count == len(indices):
229        # No points are outside
230        return array([])
231    else:
232        return indices[count:][::-1]  #return reversed
233
234
235def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
236    """Determine points inside and outside a polygon
237
238       See separate_points_by_polygon for documentation
239
240       Returns an array of points inside and an array of points outside the polygon
241    """
242
243    #if verbose: print 'Checking input to outside_polygon'
244    try:
245        points = ensure_numeric(points, Float)
246    except NameError, e:
247        raise NameError, e
248    except:
249        msg = 'Points could not be converted to Numeric array'
250        raise msg
251
252    try:
253        polygon = ensure_numeric(polygon, Float)
254    except NameError, e:
255        raise NameError, e
256    except:
257        msg = 'Polygon could not be converted to Numeric array'
258        raise msg
259
260    if len(points.shape) == 1:
261        # Only one point was passed in. Convert to array of points
262        points = reshape(points, (1,2))
263
264
265    indices, count = separate_points_by_polygon(points, polygon,
266                                                closed=closed,
267                                                verbose=verbose)
268
269    # Returns indices of points inside and indices of points outside
270    # the polygon
271
272    if count == len(indices):
273        # No points are outside
274        return indices[:count],[]
275    else:
276        return  indices[:count], indices[count:][::-1]  #return reversed
277
278
279def separate_points_by_polygon(points, polygon,
280                               closed = True, verbose = False):
281    """Determine whether points are inside or outside a polygon
282
283    Input:
284       points - Tuple of (x, y) coordinates, or list of tuples
285       polygon - list of vertices of polygon
286       closed - (optional) determine whether points on boundary should be
287       regarded as belonging to the polygon (closed = True)
288       or not (closed = False)
289
290    Outputs:
291       indices: array of same length as points with indices of points falling
292       inside the polygon listed from the beginning and indices of points
293       falling outside listed from the end.
294
295       count: count of points falling inside the polygon
296
297       The indices of points inside are obtained as indices[:count]
298       The indices of points outside are obtained as indices[count:]
299
300
301    Examples:
302       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
303
304       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
305       will return the indices [0, 2, 1] and count == 2 as only the first
306       and the last point are inside the unit square
307
308    Remarks:
309       The vertices may be listed clockwise or counterclockwise and
310       the first point may optionally be repeated.
311       Polygons do not need to be convex.
312       Polygons can have holes in them and points inside a hole is
313       regarded as being outside the polygon.
314
315    Algorithm is based on work by Darel Finley,
316    http://www.alienryderflex.com/polygon/
317
318    Uses underlying C-implementation in polygon_ext.c
319    """
320
321
322    #if verbose: print 'Checking input to separate_points_by_polygon'
323
324
325    #Input checks
326
327    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
328    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
329
330
331    try:
332        points = ensure_numeric(points, Float)
333    except NameError, e:
334        raise NameError, e
335    except:
336        msg = 'Points could not be converted to Numeric array'
337        raise msg
338
339    #if verbose: print 'Checking input to separate_points_by_polygon 2'
340    try:
341        polygon = ensure_numeric(polygon, Float)
342    except NameError, e:
343        raise NameError, e
344    except:
345        msg = 'Polygon could not be converted to Numeric array'
346        raise msg
347
348    #if verbose: print 'check'
349
350    assert len(polygon.shape) == 2,\
351       'Polygon array must be a 2d array of vertices'
352
353    assert polygon.shape[1] == 2,\
354       'Polygon array must have two columns'
355
356    assert len(points.shape) == 2,\
357       'Points array must be a 2d array'
358
359    assert points.shape[1] == 2,\
360       'Points array must have two columns'
361
362    N = polygon.shape[0] #Number of vertices in polygon
363    M = points.shape[0]  #Number of points
364
365
366    indices = zeros( M, Int )
367
368    count = _separate_points_by_polygon(points, polygon, indices,
369                                        int(closed), int(verbose))
370
371    if verbose: print 'Found %d points (out of %d) inside polygon'\
372       %(count, M)
373    return indices, count
374
375
376def polygon_area(polygon):
377    """ Determin area of arbitrary polygon
378    Reference
379    http://mathworld.wolfram.com/PolygonArea.html
380    """
381
382    n = len(polygon)
383    poly_area = 0.0
384
385    for i in range(n):
386        pti = polygon[i]
387        if i == n-1:
388            pt1 = polygon[0]
389        else:
390            pt1 = polygon[i+1]
391        xi = pti[0]
392        yi1 = pt1[1]
393        xi1 = pt1[0]
394        yi = pti[1]
395        poly_area += xi*yi1 - xi1*yi
396
397    return abs(poly_area/2)
398
399def plot_polygons(polygons_points, style=None,
400                  figname=None, label=None, verbose=False):
401
402    """ Take list of polygons and plot.
403
404    Inputs:
405
406    polygons         - list of polygons
407
408    style            - style list corresponding to each polygon
409                     - for a polygon, use 'line'
410                     - for points falling outside a polygon, use 'outside'
411
412    figname          - name to save figure to
413
414    label            - title for plot
415
416    Outputs:
417
418    - list of min and max of x and y coordinates
419    - plot of polygons
420    """
421
422    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close, title
423
424    assert type(polygons_points) == list,\
425               'input must be a list of polygons and/or points'
426
427    ion()
428    hold(True)
429
430    minx = 1e10
431    maxx = 0.0
432    miny = 1e10
433    maxy = 0.0
434
435    if label is None: label = ''
436
437    n = len(polygons_points)
438    colour = []
439    if style is None:
440        style_type = 'line'
441        style = []
442        for i in range(n):
443            style.append(style_type)
444            colour.append('b-')
445    else:
446        for s in style:
447            if s == 'line': colour.append('b-')
448            if s == 'outside': colour.append('r.')
449            if s <> 'line':
450                if s <> 'outside':
451                    colour.append('g.')
452
453    for i, item in enumerate(polygons_points):
454        x, y = poly_xy(item)
455        if min(x) < minx: minx = min(x)
456        if max(x) > maxx: maxx = max(x)
457        if min(y) < miny: miny = min(y)
458        if max(y) > maxy: maxy = max(y)
459        plot(x,y,colour[i])
460        xlabel('x')
461        ylabel('y')
462        title(label)
463
464    #raw_input('wait 1')
465    #FIXME(Ole): This makes for some strange scalings sometimes.
466    #if minx <> 0:
467    #    axis([minx*0.9,maxx*1.1,miny*0.9,maxy*1.1])
468    #else:
469    #    if miny == 0:
470    #        axis([-maxx*.01,maxx*1.1,-maxy*0.01,maxy*1.1])
471    #    else:
472    #        axis([-maxx*.01,maxx*1.1,miny*0.9,maxy*1.1])
473
474    if figname is not None:
475        savefig(figname)
476    else:
477        savefig('test_image')
478
479    close('all')
480
481    vec = [minx,maxx,miny,maxy]
482
483    return vec
484
485def poly_xy(polygon, verbose=False):
486    """ this is used within plot_polygons so need to duplicate
487        the first point so can have closed polygon in plot
488    """
489
490    #if verbose: print 'Checking input to poly_xy'
491
492    try:
493        polygon = ensure_numeric(polygon, Float)
494    except NameError, e:
495        raise NameError, e
496    except:
497        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
498        raise msg
499
500    x = polygon[:,0]
501    y = polygon[:,1]
502    x = concatenate((x, [polygon[0,0]]), axis = 0)
503    y = concatenate((y, [polygon[0,1]]), axis = 0)
504
505    return x, y
506
507#    x = []
508#    y = []
509#    n = len(poly)
510#    firstpt = poly[0]
511#    for i in range(n):
512#        thispt = poly[i]
513#        x.append(thispt[0])
514#        y.append(thispt[1])
515
516#    x.append(firstpt[0])
517#    y.append(firstpt[1])
518
519#    return x, y
520
521class Polygon_function:
522    """Create callable object f: x,y -> z, where a,y,z are vectors and
523    where f will return different values depending on whether x,y belongs
524    to specified polygons.
525
526    To instantiate:
527
528       Polygon_function(polygons)
529
530    where polygons is a list of tuples of the form
531
532      [ (P0, v0), (P1, v1), ...]
533
534      with Pi being lists of vertices defining polygons and vi either
535      constants or functions of x,y to be applied to points with the polygon.
536
537    The function takes an optional argument, default which is the value
538    (or function) to used for points not belonging to any polygon.
539    For example:
540
541       Polygon_function(polygons, default = 0.03)
542
543    If omitted the default value will be 0.0
544
545    Note: If two polygons overlap, the one last in the list takes precedence
546
547    Coordinates specified in the call are assumed to be relative to the
548    origin (georeference) e.g. used by domain.
549    By specifying the optional argument georeference,
550    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).'
561            msg += 'Got %s' %polygons
562            raise msg
563
564
565        T = regions[0]
566        try:
567            a = len(T)
568        except:
569            msg = 'Polygon_function takes a list of pairs (polygon, value).'
570            msg += 'Got %s' %polygons
571            raise msg
572
573        assert a == 2, 'Must have two component each: %s' %T
574
575
576        if geo_reference is None:
577            from anuga.coordinate_transforms.geo_reference import Geo_reference
578            geo_reference = Geo_reference()
579
580
581        self.default = default
582
583        # Make points in polygons relative to geo_reference
584        self.regions = []
585        for polygon, value in regions:
586            P = geo_reference.change_points_geo_ref(polygon)
587            self.regions.append( (P, value) )
588
589
590
591
592    def __call__(self, x, y):
593        x = array(x).astype(Float)
594        y = array(y).astype(Float)
595
596        N = len(x)
597        assert len(y) == N
598
599        points = concatenate( (reshape(x, (N, 1)),
600                               reshape(y, (N, 1))), axis=1 )
601
602        if callable(self.default):
603            z = self.default(x,y)
604        else:
605            z = ones(N, Float) * self.default
606
607        for polygon, value in self.regions:
608            indices = inside_polygon(points, polygon)
609
610            # FIXME: This needs to be vectorised
611            if callable(value):
612                for i in indices:
613                    xx = array([x[i]])
614                    yy = array([y[i]])
615                    z[i] = value(xx, yy)[0]
616            else:
617                for i in indices:
618                    z[i] = value
619
620        return z
621
622
624    """Read points assumed to form a polygon.
625       There must be exactly two numbers in each line separated by a comma.
627    """
628
629    #Get polygon
630    fid = open(filename)
632    fid.close()
633    polygon = []
634    for line in lines:
635        fields = line.split(split)
636        polygon.append( [float(fields[0]), float(fields[1])] )
637
638    return polygon
639
640
641def write_polygon(polygon, filename=None):
642    """Write polygon to csv file.
643       There will be exactly two numbers, easting and northing,
644       in each line separated by a comma.
645
647    """
648
649    fid = open(filename, 'w')
650    for point in polygon:
651        fid.write('%f, %f\n' %point)
652    fid.close()
653
654
655def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
656    """Populate given polygon with uniformly distributed points.
657
658    Input:
659       polygon - list of vertices of polygon
660       number_of_points - (optional) number of points
661       seed - seed for random number generator (default=None)
662       exclude - list of polygons (inside main polygon) from where points should be excluded
663
664    Output:
665       points - list of points inside polygon
666
667    Examples:
668       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
669       will return five randomly selected points inside the unit square
670    """
671
672    from random import uniform, seed as seed_function
673
674    seed_function(seed)
675
676    points = []
677
678    # Find outer extent of polygon
679    max_x = min_x = polygon[0][0]
680    max_y = min_y = polygon[0][1]
681    for point in polygon[1:]:
682        x = point[0]
683        if x > max_x: max_x = x
684        if x < min_x: min_x = x
685        y = point[1]
686        if y > max_y: max_y = y
687        if y < min_y: min_y = y
688
689
690    while len(points) < number_of_points:
691        x = uniform(min_x, max_x)
692        y = uniform(min_y, max_y)
693
694        append = False
695        if is_inside_polygon([x,y], polygon):
696
697            append = True
698
699            #Check exclusions
700            if exclude is not None:
701                for ex_poly in exclude:
702                    if is_inside_polygon([x,y], ex_poly):
703                        append = False
704
705
706        if append is True:
707            points.append([x,y])
708
709    return points
710
711
712def point_in_polygon(polygon, delta=1e-8):
713    """Return a point inside a given polygon which will be close to the
714    polygon edge.
715
716    Input:
717       polygon - list of vertices of polygon
718       delta - the square root of 2 * delta is the maximum distance from the
719       polygon points and the returned point.
720    Output:
721       points - a point inside polygon
722
723       searches in all diagonals and up and down (not left and right)
724    """
725    import exceptions
726    class Found(exceptions.Exception): pass
727
728    point_in = False
729    while not point_in:
730        try:
731            for poly_point in polygon: #[1:]:
732                for x_mult in range (-1,2):
733                    for y_mult in range (-1,2):
734                        x = poly_point[0]
735                        y = poly_point[1]
736                        if x == 0:
737                            x_delta = x_mult*delta
738                        else:
739                            x_delta = x+x_mult*x*delta
740
741                        if y == 0:
742                            y_delta = y_mult*delta
743                        else:
744                            y_delta = y+y_mult*y*delta
745
746                        point = [x_delta, y_delta]
747                        #print "point",point
748                        if is_inside_polygon(point, polygon, closed=False):
749                            raise Found
750        except Found:
751            point_in = True
752        else:
753            delta = delta*0.1
754    return point
755
756
757def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
758    """Calculate the approximate number of triangles inside the
759    bounding polygon and the other interior regions
760
761    Polygon areas are converted to square Kms
762
763    FIXME: Add tests for this function
764    """
765
766    from anuga.utilities.polygon import polygon_area
767
768
769    # TO DO check if any of the regions fall inside one another
770
771    print '----------------------------------------------------------------------------'
772    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
773    print '----------------------------------------------------------------------------'
774
775    no_triangles = 0.0
776    area = polygon_area(bounding_poly)
777
778    for poly, resolution in interior_regions:
779        this_area = polygon_area(poly)
780        this_triangles = this_area/resolution
781        no_triangles += this_triangles
782        area -= this_area
783
784        print 'Interior ',
785        print ('%.0f' %resolution).ljust(25),
786        print ('%.2f' %(this_area/1000000)).ljust(19),
787        print '%d' %(this_triangles)
788
789    bound_triangles = area/remainder_res
790    no_triangles += bound_triangles
791
792    print 'Bounding ',
793    print ('%.0f' %remainder_res).ljust(25),
794    print ('%.2f' %(area/1000000)).ljust(19),
795    print '%d' %(bound_triangles)
796
797    total_number_of_triangles = no_triangles/0.7
798
799    print 'Estimated total number of triangles: %d' %total_number_of_triangles
800    print 'Note: This is generally about 20% less than the final amount'
801
802    return int(total_number_of_triangles)
803
804
805def decimate_polygon(polygon, factor=10):
806    """Reduce number of points in polygon by the specified
807    factor (default=10, hence the name of the function) such that
808    the extrema in both axes are preserved.
809
810    Return reduced polygon
811    """
812
813    # FIXME(Ole): This doesn't work at present,
814    # but it isn't critical either
815
816    # Find outer extent of polygon
817    num_polygon = ensure_numeric(polygon)
818    max_x = max(num_polygon[:,0])
819    max_y = max(num_polygon[:,1])
820    min_x = min(num_polygon[:,0])
821    min_y = min(num_polygon[:,1])
822
823    # Keep only some points making sure extrema are kept
824    reduced_polygon = []
825    for i, point in enumerate(polygon):
826        x = point[0]
827        y = point[1]
828        if x in [min_x, max_x] and y in [min_y, max_y]:
829            # Keep
830            reduced_polygon.append(point)
831        else:
832            if len(reduced_polygon)*factor < i:
833                reduced_polygon.append(point)
834
835    return reduced_polygon
836
837##############################################
838#Initialise module
839
840from anuga.utilities.compile import can_use_C_extension
841if can_use_C_extension('polygon_ext.c'):
842    # Underlying C implementations can be accessed
843    from polygon_ext import _point_on_line
844    from polygon_ext import _separate_points_by_polygon
845
846else:
847    msg = 'C implementations could not be accessed by %s.\n ' %__file__
848    msg += 'Make sure compile_all.py has been run as described in '
849    msg += 'the ANUGA installation guide.'
850    raise Exception, msg
851
852
853if __name__ == "__main__":
854    pass
Note: See TracBrowser for help on using the repository browser.