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

Last change on this file since 4856 was 4856, checked in by sexton, 16 years ago

rename plot_polygons_points to original name of plot_polygons

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