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

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

Removed the less important output for fit and polygon to address ticket:238

File size: 20.6 KB
Line 
1#!/usr/bin/env python
2"""Polygon manipulations
3
4"""
5
6
7try:
8    from scipy import Float, Int, zeros, ones, array, concatenate, reshape, dot
9except:
10    #print 'Could not find scipy - using Numeric'
11    from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot
12
13
14from math import sqrt
15from anuga.utilities.numerical_tools import ensure_numeric
16from anuga.geospatial_data.geospatial_data import ensure_absolute
17
18def point_on_line(x, y, x0, y0, x1, y1):
19    """Determine whether a point is on a line segment
20
21    Input: x, y, x0, x0, x1, y1: where
22        point is given by x, y
23        line is given by (x0, y0) and (x1, y1)
24
25    """
26
27    a = array([x - x0, y - y0])
28    a_normal = array([a[1], -a[0]])
29
30    b = array([x1 - x0, y1 - y0])
31
32    if dot(a_normal, b) == 0:
33        #Point is somewhere on the infinite extension of the line
34
35        len_a = sqrt(sum(a**2))
36        len_b = sqrt(sum(b**2))
37        if dot(a, b) >= 0 and len_a <= len_b:
38           return True
39        else:
40           return False
41    else:
42      return False
43
44
45def is_inside_polygon(point, polygon, closed=True, verbose=False):
46    """Determine if one point is inside a polygon
47
48    See inside_polygon for more details
49    """
50
51    indices = inside_polygon(point, polygon, closed, verbose)
52
53    if indices.shape[0] == 1:
54        return True
55    elif indices.shape[0] == 0:
56        return False
57    else:
58        msg = 'is_inside_polygon must be invoked with one point only'
59        raise msg
60   
61
62def inside_polygon(points, polygon, closed=True, verbose=False):
63    """Determine points inside a polygon
64
65       Functions inside_polygon and outside_polygon have been defined in
66       terms af separate_by_polygon which will put all inside indices in
67       the first part of the indices array and outside indices in the last
68
69       See separate_points_by_polygon for documentation
70
71       points and polygon can be a geospatial instance,
72       a list or a numeric array
73    """
74
75    #if verbose: print 'Checking input to inside_polygon'
76
77    try:
78        points = ensure_absolute(points)
79    except NameError, e:
80        raise NameError, e
81    except:
82        # 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[0] == 1:
121        return True
122    elif indices.shape[0] == 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    Uses underlying C-implementation in polygon_ext.c
257    """
258
259
260    #if verbose: print 'Checking input to separate_points_by_polygon'
261
262
263    #Input checks
264
265    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
266    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
267
268
269    try:
270        points = ensure_numeric(points, Float)
271    except NameError, e:
272        raise NameError, e
273    except:
274        msg = 'Points could not be converted to Numeric array'
275        raise msg
276
277    #if verbose: print 'Checking input to separate_points_by_polygon 2'
278    try:
279        polygon = ensure_numeric(polygon, Float)
280    except NameError, e:
281        raise NameError, e
282    except:
283        msg = 'Polygon could not be converted to Numeric array'
284        raise msg
285
286    #if verbose: print 'check'
287
288    assert len(polygon.shape) == 2,\
289       'Polygon array must be a 2d array of vertices'
290
291    assert polygon.shape[1] == 2,\
292       'Polygon array must have two columns'
293
294    assert len(points.shape) == 2,\
295       'Points array must be a 2d array'
296
297    assert points.shape[1] == 2,\
298       'Points array must have two columns'
299
300    N = polygon.shape[0] #Number of vertices in polygon
301    M = points.shape[0]  #Number of points
302
303
304    indices = zeros( M, Int )
305
306    from polygon_ext import separate_points_by_polygon as sep_points
307    count = sep_points(points, polygon, indices,
308                       int(closed), int(verbose))
309
310    if verbose: print 'Found %d points (out of %d) inside polygon'\
311       %(count, M)
312    return indices, count
313
314
315def polygon_area(polygon):
316    """ Determin area of arbitrary polygon
317    Reference
318    http://mathworld.wolfram.com/PolygonArea.html
319    """
320   
321    n = len(polygon)
322    poly_area = 0.0
323
324    for i in range(n):
325        pti = polygon[i]
326        if i == n-1:
327            pt1 = polygon[0]
328        else:
329            pt1 = polygon[i+1]
330        xi = pti[0]
331        yi1 = pt1[1]
332        xi1 = pt1[0]
333        yi = pti[1]
334        poly_area += xi*yi1 - xi1*yi
335       
336    return abs(poly_area/2)
337
338def plot_polygons(polygons_points, style=None, 
339                         figname=None, label=None, verbose=False):
340   
341    """ Take list of polygons and plot.
342
343    Inputs:
344
345    polygons         - list of polygons
346
347    style            - style list corresponding to each polygon
348                     - for a polygon, use 'line'
349                     - for points falling outside a polygon, use 'outside'
350                       
351    figname          - name to save figure to
352
353    label            - title for plot
354
355    Outputs:
356
357    - list of min and max of x and y coordinates
358    - plot of polygons
359    """ 
360
361    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close, title
362
363    assert type(polygons_points) == list,\
364               'input must be a list of polygons and/or points'
365               
366    ion()
367    hold(True)
368
369    minx = 1e10
370    maxx = 0.0
371    miny = 1e10
372    maxy = 0.0
373
374    if label is None: label = ''
375
376    n = len(polygons_points)
377    colour = []
378    if style is None:
379        style_type = 'line' 
380        style = []
381        for i in range(n):
382            style.append(style_type)
383            colour.append('b-')
384    else:
385        for s in style:
386            if s == 'line': colour.append('b-')           
387            if s == 'outside': colour.append('r.')
388            if s <> 'line':
389                if s <> 'outside':
390                    colour.append('g.')
391           
392    for i, item in enumerate(polygons_points):
393        x, y = poly_xy(item)
394        if min(x) < minx: minx = min(x)
395        if max(x) > maxx: maxx = max(x)
396        if min(y) < miny: miny = min(y)
397        if max(y) > maxy: maxy = max(y)
398        plot(x,y,colour[i])
399        xlabel('x')
400        ylabel('y')
401        title(label)
402       
403    if minx <> 0:
404        axis([minx*0.9,maxx*1.1,miny*0.9,maxy*1.1])
405    else:
406        if miny == 0:
407            axis([-maxx*.01,maxx*1.1,-maxy*0.01,maxy*1.1])
408        else:
409            axis([-maxx*.01,maxx*1.1,miny*0.9,maxy*1.1])
410   
411    if figname is not None:
412        savefig(figname)
413    else:
414        savefig('test_image')
415
416    close('all')
417
418    vec = [minx,maxx,miny,maxy]
419
420    return vec
421
422def poly_xy(polygon, verbose=False):
423    """ this is used within plot_polygons so need to duplicate
424        the first point so can have closed polygon in plot
425    """
426
427    #if verbose: print 'Checking input to poly_xy'
428
429    try:
430        polygon = ensure_numeric(polygon, Float)
431    except NameError, e:
432        raise NameError, e
433    except:
434        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
435        raise msg
436
437    x = polygon[:,0]
438    y = polygon[:,1]
439    x = concatenate((x, [polygon[0,0]]), axis = 0)
440    y = concatenate((y, [polygon[0,1]]), axis = 0)
441   
442    return x, y
443   
444#    x = []
445#    y = []
446#    n = len(poly)
447#    firstpt = poly[0]
448#    for i in range(n):
449#        thispt = poly[i]
450#        x.append(thispt[0])
451#        y.append(thispt[1])
452
453#    x.append(firstpt[0])
454#    y.append(firstpt[1])
455   
456#    return x, y
457
458class Polygon_function:
459    """Create callable object f: x,y -> z, where a,y,z are vectors and
460    where f will return different values depending on whether x,y belongs
461    to specified polygons.
462
463    To instantiate:
464
465       Polygon_function(polygons)
466
467    where polygons is a list of tuples of the form
468
469      [ (P0, v0), (P1, v1), ...]
470
471      with Pi being lists of vertices defining polygons and vi either
472      constants or functions of x,y to be applied to points with the polygon.
473
474    The function takes an optional argument, default which is the value
475    (or function) to used for points not belonging to any polygon.
476    For example:
477
478       Polygon_function(polygons, default = 0.03)
479
480    If omitted the default value will be 0.0
481
482    Note: If two polygons overlap, the one last in the list takes precedence
483
484    Coordinates specified in the call are assumed to be relative to the origin (georeference)
485    e.g. used by domain. By specifying the optional argument georeference, all points are made relative.
486
487    FIXME: This should really work with geo_spatial point sets.
488    """
489
490    def __init__(self, regions, default = 0.0, geo_reference = None):
491
492        try:
493            len(regions)
494        except:
495            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
496            raise msg
497
498
499        T = regions[0]
500        try:
501            a = len(T)
502        except:
503            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
504            raise msg
505
506        assert a == 2, 'Must have two component each: %s' %T
507
508
509        if geo_reference is None:
510            from anuga.coordinate_transforms.geo_reference import Geo_reference
511            geo_reference = Geo_reference()
512
513
514        self.default = default
515
516        #Make points in polygons relative to geo_reference
517        self.regions = []
518        for polygon, value in regions:
519            P = geo_reference.change_points_geo_ref(polygon)
520            self.regions.append( (P, value) )
521
522
523
524
525    def __call__(self, x, y):
526        x = array(x).astype(Float)
527        y = array(y).astype(Float)
528
529        N = len(x)
530        assert len(y) == N
531
532        points = concatenate( (reshape(x, (N, 1)),
533                               reshape(y, (N, 1))), axis=1 )
534
535        if callable(self.default):
536            z = self.default(x,y)
537        else:
538            z = ones(N, Float) * self.default
539
540        for polygon, value in self.regions:
541            indices = inside_polygon(points, polygon)
542
543            #FIXME: This needs to be vectorised
544            if callable(value):
545                for i in indices:
546                    xx = array([x[i]])
547                    yy = array([y[i]])
548                    z[i] = value(xx, yy)[0]
549            else:
550                for i in indices:
551                    z[i] = value
552
553        return z
554
555
556def read_polygon(filename, split=','):
557    """Read points assumed to form a polygon.
558       There must be exactly two numbers in each line separated by a comma.
559       No header.
560    """
561
562    #Get polygon
563    fid = open(filename)
564    lines = fid.readlines()
565    fid.close()
566    polygon = []
567    for line in lines:
568        fields = line.split(split)
569        polygon.append( [float(fields[0]), float(fields[1])] )
570
571    return polygon
572
573
574def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
575    """Populate given polygon with uniformly distributed points.
576
577    Input:
578       polygon - list of vertices of polygon
579       number_of_points - (optional) number of points
580       seed - seed for random number generator (default=None)
581       exclude - list of polygons (inside main polygon) from where points should be excluded
582
583    Output:
584       points - list of points inside polygon
585
586    Examples:
587       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
588       will return five randomly selected points inside the unit square
589    """
590
591    from random import uniform, seed as seed_function
592
593    seed_function(seed)
594
595    points = []
596
597    #Find outer extent of polygon
598    max_x = min_x = polygon[0][0]
599    max_y = min_y = polygon[0][1]
600    for point in polygon[1:]:
601        x = point[0]
602        if x > max_x: max_x = x
603        if x < min_x: min_x = x
604        y = point[1]
605        if y > max_y: max_y = y
606        if y < min_y: min_y = y
607
608
609    while len(points) < number_of_points:
610        x = uniform(min_x, max_x)
611        y = uniform(min_y, max_y)
612
613        append = False
614        if is_inside_polygon([x,y], polygon):
615
616            append = True
617
618            #Check exclusions
619            if exclude is not None:
620                for ex_poly in exclude:
621                    if is_inside_polygon([x,y], ex_poly):
622                        append = False
623
624
625        if append is True:
626            points.append([x,y])
627
628    return points
629
630
631def point_in_polygon(polygon, delta=1e-8):
632    """Return a point inside a given polygon which will be close to the
633    polygon edge.
634
635    Input:
636       polygon - list of vertices of polygon
637       delta - the square root of 2 * delta is the maximum distance from the
638       polygon points and the returned point.
639    Output:
640       points - a point inside polygon
641
642       searches in all diagonals and up and down (not left and right)
643    """
644    import exceptions
645    class Found(exceptions.Exception): pass
646
647    point_in = False
648    while not point_in:
649        try:
650            for poly_point in polygon: #[1:]:
651                for x_mult in range (-1,2):
652                    for y_mult in range (-1,2):
653                        x = poly_point[0]
654                        y = poly_point[1]
655                        if x == 0:
656                            x_delta = x_mult*delta
657                        else:
658                            x_delta = x+x_mult*x*delta
659
660                        if y == 0:
661                            y_delta = y_mult*delta
662                        else:
663                            y_delta = y+y_mult*y*delta
664
665                        point = [x_delta, y_delta]
666                        #print "point",point
667                        if is_inside_polygon(point, polygon, closed=False):
668                            raise Found
669        except Found:
670            point_in = True
671        else:
672            delta = delta*0.1
673    return point
674
675
676def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
677    """Calculate the approximate number of triangles inside the
678    bounding polygon and the other interior regions
679
680    Polygon areas are converted to square Kms
681
682    FIXME: Add tests for this function
683    """
684   
685    from anuga.utilities.polygon import polygon_area
686
687
688    # TO DO check if any of the regions fall inside one another
689
690    print '----------------------------------------------------------------------------'
691    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
692    print '----------------------------------------------------------------------------'   
693       
694    no_triangles = 0.0
695    area = polygon_area(bounding_poly)
696   
697    for poly, resolution in interior_regions:
698        this_area = polygon_area(poly)
699        this_triangles = this_area/resolution
700        no_triangles += this_triangles
701        area -= this_area
702       
703        print 'Interior ',
704        print ('%.0f' %resolution).ljust(25),
705        print ('%.2f' %(this_area/1000000)).ljust(19),
706        print '%d' %(this_triangles)
707       
708    bound_triangles = area/remainder_res
709    no_triangles += bound_triangles
710
711    print 'Bounding ',
712    print ('%.0f' %remainder_res).ljust(25),
713    print ('%.2f' %(area/1000000)).ljust(19),
714    print '%d' %(bound_triangles)   
715
716    total_number_of_triangles = no_triangles/0.7
717
718    print 'Estimated total number of triangles: %d' %total_number_of_triangles
719    print 'Note: This is generally about 20% less than the final amount'   
720
721    return int(total_number_of_triangles)
722
723
724##############################################
725#Initialise module
726
727from anuga.utilities.compile import can_use_C_extension
728if can_use_C_extension('polygon_ext.c'):
729    # Underlying C implementations can be accessed
730
731    from polygon_ext import point_on_line
732else:
733    msg = 'C implementations could not be accessed by %s.\n ' %__file__
734    msg += 'Make sure compile_all.py has been run as described in '
735    msg += 'the ANUGA installation guide.'
736    raise Exception, msg
737
738
739if __name__ == "__main__":
740    pass
Note: See TracBrowser for help on using the repository browser.