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

Last change on this file since 4095 was 4095, checked in by nick, 17 years ago

add number of trigs in mesh

File size: 20.9 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        # FIXME(Ole): This message is wrong. Shouldn't it be "Absolute"
83        # rather than "Numeric"?
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        # FIXME(Ole): This message is wrong. Shouldn't it be "Absolute"
93        # rather than "Numeric"?
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    """
257
258    #Input checks
259
260
261    try:
262        points = ensure_numeric(points, Float)
263    except NameError, e:
264        raise NameError, e
265    except:
266        msg = 'Points could not be converted to Numeric array'
267        raise msg
268
269    try:
270        polygon = ensure_numeric(polygon, Float)
271    except NameError, e:
272        raise NameError, e
273    except:
274        msg = 'Polygon could not be converted to Numeric array'
275        raise msg
276
277    assert len(polygon.shape) == 2,\
278       'Polygon array must be a 2d array of vertices'
279
280    assert polygon.shape == 2,\
281       'Polygon array must have two columns'
282
283    assert len(points.shape) == 2,\
284       'Points array must be a 2d array'
285
286    assert points.shape == 2,\
287       'Points array must have two columns'
288
289    N = polygon.shape #Number of vertices in polygon
290    M = points.shape  #Number of points
291
292    px = polygon[:,0]
293    py = polygon[:,1]
294
295    #Used for an optimisation when points are far away from polygon
296    minpx = min(px); maxpx = max(px)
297    minpy = min(py); maxpy = max(py)
298
299
300    #Begin main loop
301    indices = zeros(M, Int)
302
303    inside_index = 0    #Keep track of points inside
304    outside_index = M-1 #Keep track of points outside (starting from end)
305
306    if verbose: print 'Separating %d points' %M
307    for k in range(M):
308
309        if verbose:
310            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
311
312        #for each point
313        x = points[k, 0]
314        y = points[k, 1]
315
316        inside = False
317
318        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
319            #Check polygon
320            for i in range(N):
321                j = (i+1)%N
322
323                #Check for case where point is contained in line segment
324                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
325                    if closed:
326                        inside = True
327                    else:
328                        inside = False
329                    break
330                else:
331                    #Check if truly inside polygon
332                    if py[i] < y and py[j] >= y or\
333                       py[j] < y and py[i] >= y:
334                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
335                            inside = not inside
336
337        if inside:
338            indices[inside_index] = k
339            inside_index += 1
340        else:
341            indices[outside_index] = k
342            outside_index -= 1
343
344    return indices, inside_index
345
346
347def separate_points_by_polygon_c(points, polygon,
348                                 closed = True, verbose = False):
349    """Determine whether points are inside or outside a polygon
350
351    C-wrapper
352    """
353
354
355    if verbose: print 'Checking input to separate_points_by_polygon'
356
357
358    #Input checks
359
360    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
361    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
362
363
364    try:
365        points = ensure_numeric(points, Float)
366    except NameError, e:
367        raise NameError, e
368    except:
369        msg = 'Points could not be converted to Numeric array'
370        raise msg
371
372    #if verbose: print 'Checking input to separate_points_by_polygon 2'
373    try:
374        polygon = ensure_numeric(polygon, Float)
375    except NameError, e:
376        raise NameError, e
377    except:
378        msg = 'Polygon could not be converted to Numeric array'
379        raise msg
380
381    #if verbose: print 'check'
382
383    assert len(polygon.shape) == 2,\
384       'Polygon array must be a 2d array of vertices'
385
386    assert polygon.shape == 2,\
387       'Polygon array must have two columns'
388
389    assert len(points.shape) == 2,\
390       'Points array must be a 2d array'
391
392    assert points.shape == 2,\
393       'Points array must have two columns'
394
395    N = polygon.shape #Number of vertices in polygon
396    M = points.shape  #Number of points
397
398    from polygon_ext import separate_points_by_polygon
399
400    if verbose: print 'Allocating array for indices'
401
402    indices = zeros( M, Int )
403
404    #if verbose: print 'Calling C-version of inside poly'
405
406    count = separate_points_by_polygon(points, polygon, indices,
407                                       int(closed), int(verbose))
408
409    if verbose: print 'Found %d points (out of %d) inside polygon'\
410       %(count, M)
411    return indices, count
412
413
414def polygon_area(polygon):
415    """ Determin area of arbitrary polygon
416    Reference
417    http://mathworld.wolfram.com/PolygonArea.html
418    """
419
420    n = len(polygon)
421    poly_area = 0.0
422
423    for i in range(n):
424        pti = polygon[i]
425        if i == n-1:
426            pt1 = polygon
427        else:
428            pt1 = polygon[i+1]
429        xi = pti
430        yi1 = pt1
431        xi1 = pt1
432        yi = pti
433        poly_area += xi*yi1 - xi1*yi
434
435    return abs(poly_area/2)
436
437def plot_polygons(polygons, figname, verbose=False):
438
439    """ Take list of polygons and plot.
440
441    Inputs:
442
443    polygons         - list of polygons
444
445    figname          - name to save figure to
446
447    Outputs:
448
449    - list of min and max of x and y coordinates
450    - plot of polygons
451    """
452
453    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close
454
455    assert type(polygons) == list,\
456               'input must be a list of polygons'
457
458    ion()
459    hold(True)
460
461    minx = 1e10
462    maxx = 0.0
463    miny = 1e10
464    maxy = 0.0
465
466    for polygon in polygons:
467        x, y = poly_xy(polygon)
468        if min(x) < minx: minx = min(x)
469        if max(x) > maxx: maxx = max(x)
470        if min(y) < miny: miny = min(y)
471        if max(y) > maxy: maxy = max(y)
472        plot(x,y,'r-')
473        xlabel('x')
474        ylabel('y')
475
476    savefig(figname)
477
478    close('all')
479
480    vec = [minx,maxx,miny,maxy]
481
482    return vec
483
484def poly_xy(polygon, verbose=False):
485    """ this is used within plot_polygons so need to duplicate
486        the first point so can have closed polygon in plot
487    """
488
489    if verbose: print 'Checking input to poly_xy'
490
491    try:
492        polygon = ensure_numeric(polygon, Float)
493    except NameError, e:
494        raise NameError, e
495    except:
496        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
497        raise msg
498
499    x = polygon[:,0]
500    y = polygon[:,1]
501    x = concatenate((x, [polygon[0,0]]), axis = 0)
502    y = concatenate((y, [polygon[0,1]]), axis = 0)
503
504    return x, y
505
506#    x = []
507#    y = []
508#    n = len(poly)
509#    firstpt = poly
510#    for i in range(n):
511#        thispt = poly[i]
512#        x.append(thispt)
513#        y.append(thispt)
514
515#    x.append(firstpt)
516#    y.append(firstpt)
517
518#    return x, y
519
520class Polygon_function:
521    """Create callable object f: x,y -> z, where a,y,z are vectors and
522    where f will return different values depending on whether x,y belongs
523    to specified polygons.
524
525    To instantiate:
526
527       Polygon_function(polygons)
528
529    where polygons is a list of tuples of the form
530
531      [ (P0, v0), (P1, v1), ...]
532
533      with Pi being lists of vertices defining polygons and vi either
534      constants or functions of x,y to be applied to points with the polygon.
535
536    The function takes an optional argument, default which is the value
537    (or function) to used for points not belonging to any polygon.
538    For example:
539
540       Polygon_function(polygons, default = 0.03)
541
542    If omitted the default value will be 0.0
543
544    Note: If two polygons overlap, the one last in the list takes precedence
545
546    Coordinates specified in the call are assumed to be relative to the origin (georeference)
547    e.g. used by domain. By specifying the optional argument georeference, all points are made relative.
548
549    FIXME: This should really work with geo_spatial point sets.
550    """
551
552    def __init__(self, regions, default = 0.0, geo_reference = None):
553
554        try:
555            len(regions)
556        except:
557            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
558            raise msg
559
560
561        T = regions
562        try:
563            a = len(T)
564        except:
565            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
566            raise msg
567
568        assert a == 2, 'Must have two component each: %s' %T
569
570
571        if geo_reference is None:
572            from anuga.coordinate_transforms.geo_reference import Geo_reference
573            geo_reference = Geo_reference()
574
575
576        self.default = default
577
578        #Make points in polygons relative to geo_reference
579        self.regions = []
580        for polygon, value in regions:
581            P = geo_reference.change_points_geo_ref(polygon)
582            self.regions.append( (P, value) )
583
584
585
586
587    def __call__(self, x, y):
588        x = array(x).astype(Float)
589        y = array(y).astype(Float)
590
591        N = len(x)
592        assert len(y) == N
593
594        points = concatenate( (reshape(x, (N, 1)),
595                               reshape(y, (N, 1))), axis=1 )
596
597        if callable(self.default):
598            z = self.default(x,y)
599        else:
600            z = ones(N, Float) * self.default
601
602        for polygon, value in self.regions:
603            indices = inside_polygon(points, polygon)
604
605            #FIXME: This needs to be vectorised
606            if callable(value):
607                for i in indices:
608                    xx = array([x[i]])
609                    yy = array([y[i]])
610                    z[i] = value(xx, yy)
611            else:
612                for i in indices:
613                    z[i] = value
614
615        return z
616
618    """Read points assumed to form a polygon.
619       There must be exactly two numbers in each line separated by a comma.
621    """
622
623    #Get polygon
624    fid = open(filename)
626    fid.close()
627    polygon = []
628    for line in lines:
629        fields = line.split(split)
630        polygon.append( [float(fields), float(fields)] )
631
632    return polygon
633
634def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
635    """Populate given polygon with uniformly distributed points.
636
637    Input:
638       polygon - list of vertices of polygon
639       number_of_points - (optional) number of points
640       seed - seed for random number generator (default=None)
641       exclude - list of polygons (inside main polygon) from where points should be excluded
642
643    Output:
644       points - list of points inside polygon
645
646    Examples:
647       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
648       will return five randomly selected points inside the unit square
649    """
650
651    from random import uniform, seed as seed_function
652
653    seed_function(seed)
654
655    points = []
656
657    #Find outer extent of polygon
658    max_x = min_x = polygon
659    max_y = min_y = polygon
660    for point in polygon[1:]:
661        x = point
662        if x > max_x: max_x = x
663        if x < min_x: min_x = x
664        y = point
665        if y > max_y: max_y = y
666        if y < min_y: min_y = y
667
668
669    while len(points) < number_of_points:
670        x = uniform(min_x, max_x)
671        y = uniform(min_y, max_y)
672
673        append = False
674        if is_inside_polygon([x,y], polygon):
675
676            append = True
677
678            #Check exclusions
679            if exclude is not None:
680                for ex_poly in exclude:
681                    if is_inside_polygon([x,y], ex_poly):
682                        append = False
683
684
685        if append is True:
686            points.append([x,y])
687
688    return points
689
690def point_in_polygon(polygon, delta=1e-8):
691    """Return a point inside a given polygon which will be close to the
692    polygon edge.
693
694    Input:
695       polygon - list of vertices of polygon
696       delta - the square root of 2 * delta is the maximum distance from the
697       polygon points and the returned point.
698    Output:
699       points - a point inside polygon
700
701       searches in all diagonals and up and down (not left and right)
702    """
703    import exceptions
704    class Found(exceptions.Exception): pass
705
706    point_in = False
707    while not point_in:
708        try:
709            for poly_point in polygon: #[1:]:
710                for x_mult in range (-1,2):
711                    for y_mult in range (-1,2):
712                        x = poly_point
713                        y = poly_point
714                        if x == 0:
715                            x_delta = x_mult*delta
716                        else:
717                            x_delta = x+x_mult*x*delta
718
719                        if y == 0:
720                            y_delta = y_mult*delta
721                        else:
722                            y_delta = y+y_mult*y*delta
723
724                        point = [x_delta, y_delta]
725                        #print "point",point
726                        if is_inside_polygon(point, polygon, closed=False):
727                            raise Found
728        except Found:
729            point_in = True
730        else:
731            delta = delta*0.1
732    return point
733
734def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
735    """Calcalutes the approximate number of triangles inside the bounding polygon
736    and the other interior regions
737    FIXME: Add tests for this function
738    """
739    import polygon_area
740
741    # TO DO check if any of the regions fall inside one another
742    no_triangles = 0.0
743    area = polygon_area(bounding_poly)
744    for i,j in interior_regions:
745        this_area = polygon_area(i)
746        no_triangles += this_area/j
747        area -= this_area
748        #convert to square Kms
749        print 'area of ',j, this_area/1000000., area/1000000.
750    no_triangles += area/remainder_res
751    return int(no_triangles/0.7)
752
753
754
755##############################################
756#Initialise module
757
758from anuga.utilities.compile import can_use_C_extension
759if can_use_C_extension('polygon_ext.c'):
760    from polygon_ext import point_on_line
761    separate_points_by_polygon = separate_points_by_polygon_c
762
763
764if __name__ == "__main__":
765    pass
Note: See TracBrowser for help on using the repository browser.