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

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

Formatting

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