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

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

Decimated polygon test file to remove potentially licensed data.
Fiddled with code.

File size: 21.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
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    if minx <> 0:
405        axis([minx*0.9,maxx*1.1,miny*0.9,maxy*1.1])
406    else:
407        if miny == 0:
408            axis([-maxx*.01,maxx*1.1,-maxy*0.01,maxy*1.1])
409        else:
410            axis([-maxx*.01,maxx*1.1,miny*0.9,maxy*1.1])
411   
412    if figname is not None:
413        savefig(figname)
414    else:
415        savefig('test_image')
416
417    close('all')
418
419    vec = [minx,maxx,miny,maxy]
420
421    return vec
422
423def poly_xy(polygon, verbose=False):
424    """ this is used within plot_polygons so need to duplicate
425        the first point so can have closed polygon in plot
426    """
427
428    #if verbose: print 'Checking input to poly_xy'
429
430    try:
431        polygon = ensure_numeric(polygon, Float)
432    except NameError, e:
433        raise NameError, e
434    except:
435        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
436        raise msg
437
438    x = polygon[:,0]
439    y = polygon[:,1]
440    x = concatenate((x, [polygon[0,0]]), axis = 0)
441    y = concatenate((y, [polygon[0,1]]), axis = 0)
442   
443    return x, y
444   
445#    x = []
446#    y = []
447#    n = len(poly)
448#    firstpt = poly[0]
449#    for i in range(n):
450#        thispt = poly[i]
451#        x.append(thispt[0])
452#        y.append(thispt[1])
453
454#    x.append(firstpt[0])
455#    y.append(firstpt[1])
456   
457#    return x, y
458
459class Polygon_function:
460    """Create callable object f: x,y -> z, where a,y,z are vectors and
461    where f will return different values depending on whether x,y belongs
462    to specified polygons.
463
464    To instantiate:
465
466       Polygon_function(polygons)
467
468    where polygons is a list of tuples of the form
469
470      [ (P0, v0), (P1, v1), ...]
471
472      with Pi being lists of vertices defining polygons and vi either
473      constants or functions of x,y to be applied to points with the polygon.
474
475    The function takes an optional argument, default which is the value
476    (or function) to used for points not belonging to any polygon.
477    For example:
478
479       Polygon_function(polygons, default = 0.03)
480
481    If omitted the default value will be 0.0
482
483    Note: If two polygons overlap, the one last in the list takes precedence
484
485    Coordinates specified in the call are assumed to be relative to the origin (georeference)
486    e.g. used by domain. By specifying the optional argument georeference, all points are made relative.
487
488    FIXME: This should really work with geo_spatial point sets.
489    """
490
491    def __init__(self, regions, default = 0.0, geo_reference = None):
492
493        try:
494            len(regions)
495        except:
496            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
497            raise msg
498
499
500        T = regions[0]
501        try:
502            a = len(T)
503        except:
504            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
505            raise msg
506
507        assert a == 2, 'Must have two component each: %s' %T
508
509
510        if geo_reference is None:
511            from anuga.coordinate_transforms.geo_reference import Geo_reference
512            geo_reference = Geo_reference()
513
514
515        self.default = default
516
517        #Make points in polygons relative to geo_reference
518        self.regions = []
519        for polygon, value in regions:
520            P = geo_reference.change_points_geo_ref(polygon)
521            self.regions.append( (P, value) )
522
523
524
525
526    def __call__(self, x, y):
527        x = array(x).astype(Float)
528        y = array(y).astype(Float)
529
530        N = len(x)
531        assert len(y) == N
532
533        points = concatenate( (reshape(x, (N, 1)),
534                               reshape(y, (N, 1))), axis=1 )
535
536        if callable(self.default):
537            z = self.default(x,y)
538        else:
539            z = ones(N, Float) * self.default
540
541        for polygon, value in self.regions:
542            indices = inside_polygon(points, polygon)
543
544            #FIXME: This needs to be vectorised
545            if callable(value):
546                for i in indices:
547                    xx = array([x[i]])
548                    yy = array([y[i]])
549                    z[i] = value(xx, yy)[0]
550            else:
551                for i in indices:
552                    z[i] = value
553
554        return z
555
556
557def read_polygon(filename, split=','):
558    """Read points assumed to form a polygon.
559       There must be exactly two numbers in each line separated by a comma.
560       No header.
561    """
562
563    #Get polygon
564    fid = open(filename)
565    lines = fid.readlines()
566    fid.close()
567    polygon = []
568    for line in lines:
569        fields = line.split(split)
570        polygon.append( [float(fields[0]), float(fields[1])] )
571
572    return polygon
573
574
575def write_polygon(polygon, filename=None):
576    """Write polygon to csv file.
577       There will be exactly two numbers, easting and northing,
578       in each line separated by a comma.
579       
580       No header.   
581    """
582
583    fid = open(filename, 'w')
584    for point in polygon:
585        fid.write('%f, %f\n' %point)
586    fid.close()
587   
588
589def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
590    """Populate given polygon with uniformly distributed points.
591
592    Input:
593       polygon - list of vertices of polygon
594       number_of_points - (optional) number of points
595       seed - seed for random number generator (default=None)
596       exclude - list of polygons (inside main polygon) from where points should be excluded
597
598    Output:
599       points - list of points inside polygon
600
601    Examples:
602       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
603       will return five randomly selected points inside the unit square
604    """
605
606    from random import uniform, seed as seed_function
607
608    seed_function(seed)
609
610    points = []
611
612    # Find outer extent of polygon
613    max_x = min_x = polygon[0][0]
614    max_y = min_y = polygon[0][1]
615    for point in polygon[1:]:
616        x = point[0]
617        if x > max_x: max_x = x
618        if x < min_x: min_x = x
619        y = point[1]
620        if y > max_y: max_y = y
621        if y < min_y: min_y = y
622
623
624    while len(points) < number_of_points:
625        x = uniform(min_x, max_x)
626        y = uniform(min_y, max_y)
627
628        append = False
629        if is_inside_polygon([x,y], polygon):
630
631            append = True
632
633            #Check exclusions
634            if exclude is not None:
635                for ex_poly in exclude:
636                    if is_inside_polygon([x,y], ex_poly):
637                        append = False
638
639
640        if append is True:
641            points.append([x,y])
642
643    return points
644
645
646def point_in_polygon(polygon, delta=1e-8):
647    """Return a point inside a given polygon which will be close to the
648    polygon edge.
649
650    Input:
651       polygon - list of vertices of polygon
652       delta - the square root of 2 * delta is the maximum distance from the
653       polygon points and the returned point.
654    Output:
655       points - a point inside polygon
656
657       searches in all diagonals and up and down (not left and right)
658    """
659    import exceptions
660    class Found(exceptions.Exception): pass
661
662    point_in = False
663    while not point_in:
664        try:
665            for poly_point in polygon: #[1:]:
666                for x_mult in range (-1,2):
667                    for y_mult in range (-1,2):
668                        x = poly_point[0]
669                        y = poly_point[1]
670                        if x == 0:
671                            x_delta = x_mult*delta
672                        else:
673                            x_delta = x+x_mult*x*delta
674
675                        if y == 0:
676                            y_delta = y_mult*delta
677                        else:
678                            y_delta = y+y_mult*y*delta
679
680                        point = [x_delta, y_delta]
681                        #print "point",point
682                        if is_inside_polygon(point, polygon, closed=False):
683                            raise Found
684        except Found:
685            point_in = True
686        else:
687            delta = delta*0.1
688    return point
689
690
691def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
692    """Calculate the approximate number of triangles inside the
693    bounding polygon and the other interior regions
694
695    Polygon areas are converted to square Kms
696
697    FIXME: Add tests for this function
698    """
699   
700    from anuga.utilities.polygon import polygon_area
701
702
703    # TO DO check if any of the regions fall inside one another
704
705    print '----------------------------------------------------------------------------'
706    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
707    print '----------------------------------------------------------------------------'   
708       
709    no_triangles = 0.0
710    area = polygon_area(bounding_poly)
711   
712    for poly, resolution in interior_regions:
713        this_area = polygon_area(poly)
714        this_triangles = this_area/resolution
715        no_triangles += this_triangles
716        area -= this_area
717       
718        print 'Interior ',
719        print ('%.0f' %resolution).ljust(25),
720        print ('%.2f' %(this_area/1000000)).ljust(19),
721        print '%d' %(this_triangles)
722       
723    bound_triangles = area/remainder_res
724    no_triangles += bound_triangles
725
726    print 'Bounding ',
727    print ('%.0f' %remainder_res).ljust(25),
728    print ('%.2f' %(area/1000000)).ljust(19),
729    print '%d' %(bound_triangles)   
730
731    total_number_of_triangles = no_triangles/0.7
732
733    print 'Estimated total number of triangles: %d' %total_number_of_triangles
734    print 'Note: This is generally about 20% less than the final amount'   
735
736    return int(total_number_of_triangles)
737
738
739def decimate_polygon(polygon, factor=10):
740    """Reduce number of points in polygon by the specified
741    factor (default=10, hence the name of the function) such that
742    the extrema in both axes are preserved.
743
744    Return reduced polygon
745    """
746
747    # FIXME(Ole): This doesn't work at present,
748    # but it isn't critical either
749
750    # Find outer extent of polygon
751    num_polygon = ensure_numeric(polygon)
752    max_x = max(num_polygon[:,0])
753    max_y = max(num_polygon[:,1])
754    min_x = min(num_polygon[:,0])
755    min_y = min(num_polygon[:,1])       
756
757    # Keep only some points making sure extrema are kept
758    reduced_polygon = []   
759    for i, point in enumerate(polygon):
760        x = point[0]
761        y = point[1]       
762        if x in [min_x, max_x] and y in [min_y, max_y]:
763            # Keep
764            reduced_polygon.append(point)
765        else:
766            if len(reduced_polygon)*factor < i:
767                reduced_polygon.append(point)               
768
769    return reduced_polygon
770
771##############################################
772#Initialise module
773
774from anuga.utilities.compile import can_use_C_extension
775if can_use_C_extension('polygon_ext.c'):
776    # Underlying C implementations can be accessed
777
778    from polygon_ext import point_on_line
779else:
780    msg = 'C implementations could not be accessed by %s.\n ' %__file__
781    msg += 'Make sure compile_all.py has been run as described in '
782    msg += 'the ANUGA installation guide.'
783    raise Exception, msg
784
785
786if __name__ == "__main__":
787    pass
Note: See TracBrowser for help on using the repository browser.