source: anuga_core/source_numpy_conversion/anuga/utilities/polygon.py @ 5948

Last change on this file since 5948 was 5948, checked in by rwilson, 15 years ago

More NumPy? changes.

File size: 30.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
12import numpy
13
14from math import sqrt
15from anuga.utilities.numerical_tools import ensure_numeric
16from anuga.geospatial_data.geospatial_data import ensure_absolute
17from anuga.config import Float, Int
18
19
20def point_on_line(point, line, rtol=1.0e-5, atol=1.0e-8):
21    """Determine whether a point is on a line segment
22
23    Input:
24        point is given by [x, y]
25        line is given by [x0, y0], [x1, y1]] or
26        the equivalent 2x2 Numeric array with each row corresponding to a point.
27
28    Output:
29
30    Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding.
31    """
32
33    point = ensure_numeric(point)
34    line = ensure_numeric(line)
35
36    res = _point_on_line(point[0], point[1],
37                         line[0,0], line[0,1],
38                         line[1,0], line[1,1],
39                         rtol, atol)
40   
41    return bool(res)
42
43
44######
45# Result functions used in intersection() below for collinear lines.
46# (p0,p1) defines line 0, (p2,p3) defines line 1.
47######
48
49# result functions for possible states
50def lines_dont_coincide(p0,p1,p2,p3):               return (3, None)
51def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2, numpy.array([p0,p1]))
52def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2, numpy.array([p2,p3]))
53def lines_overlap_same_direction(p0,p1,p2,p3):      return (2, numpy.array([p0,p3]))
54def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2, numpy.array([p2,p1]))
55def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2, numpy.array([p0,p2]))
56def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, numpy.array([p3,p1]))
57
58# this function called when an impossible state is found
59def lines_error(p1, p2, p3, p4): raise RuntimeError, "INTERNAL ERROR"
60
61#                     0s1    0e1    1s0    1e0   # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0
62collinear_result = { (False, False, False, False): lines_dont_coincide,
63                     (False, False, False, True ): lines_error,
64                     (False, False, True,  False): lines_error,
65                     (False, False, True,  True ): lines_1_fully_included_in_0,
66                     (False, True,  False, False): lines_error,
67                     (False, True,  False, True ): lines_overlap_opposite_direction2,
68                     (False, True,  True,  False): lines_overlap_same_direction2,
69                     (False, True,  True,  True ): lines_1_fully_included_in_0,
70                     (True,  False, False, False): lines_error,
71                     (True,  False, False, True ): lines_overlap_same_direction,
72                     (True,  False, True,  False): lines_overlap_opposite_direction,
73                     (True,  False, True,  True ): lines_1_fully_included_in_0,
74                     (True,  True,  False, False): lines_0_fully_included_in_1,
75                     (True,  True,  False, True ): lines_0_fully_included_in_1,
76                     (True,  True,  True,  False): lines_0_fully_included_in_1,
77                     (True,  True,  True,  True ): lines_0_fully_included_in_1
78                   }
79
80def intersection(line0, line1, rtol=1.0e-5, atol=1.0e-8):
81    """Returns intersecting point between two line segments or None
82    (if parallel or no intersection is found).
83
84    However, if parallel lines coincide partly (i.e. shara a common segment,
85    the line segment where lines coincide is returned
86   
87
88    Inputs:
89        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
90                      A line can also be a 2x2 numpy array with each row
91                      corresponding to a point.
92
93
94    Output:
95        status, value
96
97        where status and value is interpreted as follows
98       
99        status == 0: no intersection, value set to None.
100        status == 1: intersection point found and returned in value as [x,y].
101        status == 2: Collinear overlapping lines found. Value takes the form [[x0,y0], [x1,y1]].
102        status == 3: Collinear non-overlapping lines. Value set to None.
103        status == 4: Lines are parallel with a fixed distance apart. Value set to None.
104   
105    """
106
107    # FIXME (Ole): Write this in C
108
109    line0 = ensure_numeric(line0, Float)
110    line1 = ensure_numeric(line1, Float)   
111
112    x0 = line0[0,0]; y0 = line0[0,1]
113    x1 = line0[1,0]; y1 = line0[1,1]
114
115    x2 = line1[0,0]; y2 = line1[0,1]
116    x3 = line1[1,0]; y3 = line1[1,1]
117
118    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
119    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
120    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
121       
122    if numpy.allclose(denom, 0.0, rtol=rtol, atol=atol):
123        # Lines are parallel - check if they are collinear
124        if numpy.allclose([u0, u1], 0.0, rtol=rtol, atol=atol):
125            # We now know that the lines are collinear
126            state_tuple = (point_on_line([x0, y0], line1, rtol=rtol, atol=atol),
127                           point_on_line([x1, y1], line1, rtol=rtol, atol=atol),
128                           point_on_line([x2, y2], line0, rtol=rtol, atol=atol),
129                           point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
130
131            return collinear_result[state_tuple]([x0,y0],[x1,y1],[x2,y2],[x3,y3])
132        else:
133            # Lines are parallel but aren't collinear
134            return 4, None #FIXME (Ole): Add distance here instead of None
135    else:
136        # Lines are not parallel, check if they intersect
137        u0 = u0/denom
138        u1 = u1/denom       
139
140        x = x0 + u0*(x1-x0)
141        y = y0 + u0*(y1-y0)
142
143        # Sanity check - can be removed to speed up if needed
144        assert numpy.allclose(x, x2 + u1*(x3-x2), rtol=rtol, atol=atol)
145        assert numpy.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol)       
146
147        # Check if point found lies within given line segments
148        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 
149            # We have intersection
150            return 1, numpy.array([x, y])
151        else:
152            # No intersection
153            return 0, None
154
155
156def NEW_C_intersection(line0, line1):
157    #FIXME(Ole): To write in C
158    """Returns intersecting point between two line segments or None
159    (if parallel or no intersection is found).
160
161    However, if parallel lines coincide partly (i.e. shara a common segment,
162    the line segment where lines coincide is returned
163   
164
165    Inputs:
166        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
167                      A line can also be a 2x2 numeric array with each row
168                      corresponding to a point.
169
170
171    Output:
172        status, value
173
174        where status is interpreted as follows
175       
176        status == 0: no intersection with value set to None
177        status == 1: One intersection point found and returned in value as [x,y]
178        status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]]
179        status == 3: Lines would coincide but only if extended. Value set to None
180        status == 4: Lines are parallel with a fixed distance apart. Value set to None.
181   
182    """
183
184
185    line0 = ensure_numeric(line0, Float)
186    line1 = ensure_numeric(line1, Float)   
187
188    status, value = _intersection(line0[0,0], line0[0,1],
189                                  line0[1,0], line0[1,1],
190                                  line1[0,0], line1[0,1],
191                                  line1[1,0], line1[1,1])
192
193    return status, value
194
195
196
197
198def is_inside_polygon(point, polygon, closed=True, verbose=False):
199    """Determine if one point is inside a polygon
200
201    See inside_polygon for more details
202    """
203
204    indices = inside_polygon(point, polygon, closed, verbose)
205
206    if indices.shape[0] == 1:
207        return True
208    elif indices.shape[0] == 0:
209        return False
210    else:
211        msg = 'is_inside_polygon must be invoked with one point only'
212        raise msg
213   
214
215def inside_polygon(points, polygon, closed=True, verbose=False):
216    """Determine points inside a polygon
217
218       Functions inside_polygon and outside_polygon have been defined in
219       terms af separate_by_polygon which will put all inside indices in
220       the first part of the indices array and outside indices in the last
221
222       See separate_points_by_polygon for documentation
223
224       points and polygon can be a geospatial instance,
225       a list or a numeric array
226    """
227
228    #if verbose: print 'Checking input to inside_polygon'
229
230    try:
231        points = ensure_absolute(points)
232    except NameError, e:
233        raise NameError, e
234    except:
235        # If this fails it is going to be because the points can't be
236        # converted to a numeric array.
237        msg = 'Points could not be converted to Numeric array' 
238        raise TypeError, msg
239
240    try:
241        polygon = ensure_absolute(polygon)
242    except NameError, e:
243        raise NameError, e
244    except:
245        # If this fails it is going to be because the points can't be
246        # converted to a numeric array.
247        msg = 'Polygon %s could not be converted to Numeric array' % (str(polygon))
248        raise TypeError, msg
249
250    if len(points.shape) == 1:
251        # Only one point was passed in. Convert to array of points
252        points = numpy.reshape(points, (1,2))
253
254    indices, count = separate_points_by_polygon(points, polygon,
255                                                closed=closed,
256                                                verbose=verbose)
257
258    # Return indices of points inside polygon
259    return indices[:count]
260
261
262
263def is_outside_polygon(point, polygon, closed=True, verbose=False,
264                       points_geo_ref=None, polygon_geo_ref=None):
265    """Determine if one point is outside a polygon
266
267    See outside_polygon for more details
268    """
269
270    indices = outside_polygon(point, polygon, closed, verbose)
271                              #points_geo_ref, polygon_geo_ref)
272
273    if indices.shape[0] == 1:
274        return True
275    elif indices.shape[0] == 0:
276        return False
277    else:
278        msg = 'is_outside_polygon must be invoked with one point only'
279        raise msg
280   
281
282def outside_polygon(points, polygon, closed = True, verbose = False):
283    """Determine points outside a polygon
284
285       Functions inside_polygon and outside_polygon have been defined in
286       terms af separate_by_polygon which will put all inside indices in
287       the first part of the indices array and outside indices in the last
288
289       See separate_points_by_polygon for documentation
290    """
291
292    #if verbose: print 'Checking input to outside_polygon'
293    try:
294        points = ensure_numeric(points, Float)
295    except NameError, e:
296        raise NameError, e
297    except:
298        msg = 'Points could not be converted to Numeric array'
299        raise msg
300
301    try:
302        polygon = ensure_numeric(polygon, Float)
303    except NameError, e:
304        raise NameError, e
305    except:
306        msg = 'Polygon could not be converted to Numeric array'
307        raise msg
308
309
310    if len(points.shape) == 1:
311        # Only one point was passed in. Convert to array of points
312        points = numpy.reshape(points, (1,2))
313
314    indices, count = separate_points_by_polygon(points, polygon,
315                                                closed=closed,
316                                                verbose=verbose)
317
318    # Return indices of points outside polygon
319    if count == len(indices):
320        # No points are outside
321        return numpy.array([])
322    else:
323        return indices[count:][::-1]  #return reversed
324       
325
326def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
327    """Determine points inside and outside a polygon
328
329       See separate_points_by_polygon for documentation
330
331       Returns an array of points inside and an array of points outside the polygon
332    """
333
334    #if verbose: print 'Checking input to outside_polygon'
335    try:
336        points = ensure_numeric(points, Float)
337    except NameError, e:
338        raise NameError, e
339    except:
340        msg = 'Points could not be converted to Numeric array'
341        raise msg
342
343    try:
344        polygon = ensure_numeric(polygon, Float)
345    except NameError, e:
346        raise NameError, e
347    except:
348        msg = 'Polygon could not be converted to Numeric array'
349        raise msg
350
351    if len(points.shape) == 1:
352        # Only one point was passed in. Convert to array of points
353        points = numpy.reshape(points, (1,2))
354
355
356    indices, count = separate_points_by_polygon(points, polygon,
357                                                closed=closed,
358                                                verbose=verbose)
359   
360    # Returns indices of points inside and indices of points outside
361    # the polygon
362
363    if count == len(indices):
364        # No points are outside
365        return indices[:count],[]
366    else:
367        return  indices[:count], indices[count:][::-1]  #return reversed
368
369
370def separate_points_by_polygon(points, polygon,
371                               closed = True, verbose = False):
372    """Determine whether points are inside or outside a polygon
373
374    Input:
375       points - Tuple of (x, y) coordinates, or list of tuples
376       polygon - list of vertices of polygon
377       closed - (optional) determine whether points on boundary should be
378       regarded as belonging to the polygon (closed = True)
379       or not (closed = False)
380
381    Outputs:
382       indices: array of same length as points with indices of points falling
383       inside the polygon listed from the beginning and indices of points
384       falling outside listed from the end.
385
386       count: count of points falling inside the polygon
387
388       The indices of points inside are obtained as indices[:count]
389       The indices of points outside are obtained as indices[count:]
390
391
392    Examples:
393       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
394
395       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
396       will return the indices [0, 2, 1] and count == 2 as only the first
397       and the last point are inside the unit square
398
399    Remarks:
400       The vertices may be listed clockwise or counterclockwise and
401       the first point may optionally be repeated.
402       Polygons do not need to be convex.
403       Polygons can have holes in them and points inside a hole is
404       regarded as being outside the polygon.
405
406    Algorithm is based on work by Darel Finley,
407    http://www.alienryderflex.com/polygon/
408
409    Uses underlying C-implementation in polygon_ext.c
410    """
411
412
413    #if verbose: print 'Checking input to separate_points_by_polygon'
414
415
416    #Input checks
417
418    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
419    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
420
421##    print 'Before: points=%s, flags=%s' % (type(points), str(points.flags))
422    try:
423##        points = numpy.ascontiguousarray(ensure_numeric(points, Float))
424        points = ensure_numeric(points, Float)
425    except NameError, e:
426        raise NameError, e
427    except:
428        msg = 'Points could not be converted to Numeric array'
429        raise TypeError, msg
430##    print 'After: points=%s, flags=%s' % (type(points), str(points.flags))
431
432    #if verbose: print 'Checking input to separate_points_by_polygon 2'
433    try:
434##        polygon = numpy.ascontiguousarray(ensure_numeric(polygon, Float))
435        polygon = ensure_numeric(polygon, Float)
436    except NameError, e:
437        raise NameError, e
438    except:
439        msg = 'Polygon could not be converted to Numeric array'
440        raise TypeError, msg
441
442    msg = 'Polygon array must be a 2d array of vertices'
443    assert len(polygon.shape) == 2, msg
444
445    msg = 'Polygon array must have two columns' 
446    assert polygon.shape[1] == 2, msg
447
448
449    msg = 'Points array must be 1 or 2 dimensional.'
450    msg += ' I got %d dimensions' %len(points.shape)
451    assert 0 < len(points.shape) < 3, msg
452
453
454    if len(points.shape) == 1:
455        # Only one point was passed in.
456        # Convert to array of points
457        points = numpy.reshape(points, (1,2))
458
459   
460    msg = 'Point array must have two columns (x,y), '
461    msg += 'I got points.shape[1] == %d' % points.shape[1]
462    assert points.shape[1] == 2, msg
463
464       
465    msg = 'Points array must be a 2d array. I got %s' %str(points[:30])
466    assert len(points.shape) == 2, msg
467
468    msg = 'Points array must have two columns'
469    assert points.shape[1] == 2, msg
470
471
472    N = polygon.shape[0] #Number of vertices in polygon
473    M = points.shape[0]  #Number of points
474
475
476    indices = numpy.zeros( M, Int )
477
478    count = _separate_points_by_polygon(points, polygon, indices,
479                                        int(closed), int(verbose))
480
481    if verbose: print 'Found %d points (out of %d) inside polygon'\
482       %(count, M)
483    return indices, count
484
485
486def polygon_area(polygon):
487    """ Determin area of arbitrary polygon
488    Reference
489    http://mathworld.wolfram.com/PolygonArea.html
490    """
491   
492    n = len(polygon)
493    poly_area = 0.0
494
495    for i in range(n):
496        pti = polygon[i]
497        if i == n-1:
498            pt1 = polygon[0]
499        else:
500            pt1 = polygon[i+1]
501        xi = pti[0]
502        yi1 = pt1[1]
503        xi1 = pt1[0]
504        yi = pti[1]
505        poly_area += xi*yi1 - xi1*yi
506       
507    return abs(poly_area/2)
508
509def plot_polygons(polygons_points, style=None, 
510                  figname=None, label=None, verbose=False):
511   
512    """ Take list of polygons and plot.
513
514    Inputs:
515
516    polygons         - list of polygons
517
518    style            - style list corresponding to each polygon
519                     - for a polygon, use 'line'
520                     - for points falling outside a polygon, use 'outside'
521                       
522    figname          - name to save figure to
523
524    label            - title for plot
525
526    Outputs:
527
528    - list of min and max of x and y coordinates
529    - plot of polygons
530    """ 
531
532    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close, title
533
534    assert type(polygons_points) == list,\
535               'input must be a list of polygons and/or points'
536               
537    ion()
538    hold(True)
539
540    minx = 1e10
541    maxx = 0.0
542    miny = 1e10
543    maxy = 0.0
544
545    if label is None: label = ''
546
547    n = len(polygons_points)
548    colour = []
549    if style is None:
550        style_type = 'line' 
551        style = []
552        for i in range(n):
553            style.append(style_type)
554            colour.append('b-')
555    else:
556        for s in style:
557            if s == 'line': colour.append('b-')           
558            if s == 'outside': colour.append('r.')
559            if s <> 'line':
560                if s <> 'outside':
561                    colour.append('g.')
562           
563    for i, item in enumerate(polygons_points):
564        x, y = poly_xy(item)
565        if min(x) < minx: minx = min(x)
566        if max(x) > maxx: maxx = max(x)
567        if min(y) < miny: miny = min(y)
568        if max(y) > maxy: maxy = max(y)
569        plot(x,y,colour[i])
570        xlabel('x')
571        ylabel('y')
572        title(label)
573
574    #raw_input('wait 1')
575    #FIXME(Ole): This makes for some strange scalings sometimes.
576    #if minx <> 0:
577    #    axis([minx*0.9,maxx*1.1,miny*0.9,maxy*1.1])
578    #else:
579    #    if miny == 0:
580    #        axis([-maxx*.01,maxx*1.1,-maxy*0.01,maxy*1.1])
581    #    else:
582    #        axis([-maxx*.01,maxx*1.1,miny*0.9,maxy*1.1])
583
584    if figname is not None:
585        savefig(figname)
586    else:
587        savefig('test_image')
588
589    close('all')
590
591    vec = [minx,maxx,miny,maxy]
592
593    return vec
594
595def poly_xy(polygon, verbose=False):
596    """ this is used within plot_polygons so need to duplicate
597        the first point so can have closed polygon in plot
598    """
599
600    #if verbose: print 'Checking input to poly_xy'
601
602    try:
603        polygon = ensure_numeric(polygon, Float)
604    except NameError, e:
605        raise NameError, e
606    except:
607        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
608        raise msg
609
610    x = polygon[:,0]
611    y = polygon[:,1]
612    x = numpy.concatenate((x, [polygon[0,0]]), axis = 0)
613    y = numpy.concatenate((y, [polygon[0,1]]), axis = 0)
614   
615    return x, y
616   
617#    x = []
618#    y = []
619#    n = len(poly)
620#    firstpt = poly[0]
621#    for i in range(n):
622#        thispt = poly[i]
623#        x.append(thispt[0])
624#        y.append(thispt[1])
625
626#    x.append(firstpt[0])
627#    y.append(firstpt[1])
628   
629#    return x, y
630
631class Polygon_function:
632    """Create callable object f: x,y -> z, where a,y,z are vectors and
633    where f will return different values depending on whether x,y belongs
634    to specified polygons.
635
636    To instantiate:
637
638       Polygon_function(polygons)
639
640    where polygons is a list of tuples of the form
641
642      [ (P0, v0), (P1, v1), ...]
643
644      with Pi being lists of vertices defining polygons and vi either
645      constants or functions of x,y to be applied to points with the polygon.
646
647    The function takes an optional argument, default which is the value
648    (or function) to used for points not belonging to any polygon.
649    For example:
650
651       Polygon_function(polygons, default = 0.03)
652
653    If omitted the default value will be 0.0
654
655    Note: If two polygons overlap, the one last in the list takes precedence
656
657    Coordinates specified in the call are assumed to be relative to the
658    origin (georeference) e.g. used by domain.
659    By specifying the optional argument georeference,
660    all points are made relative.
661
662    FIXME: This should really work with geo_spatial point sets.
663    """
664
665    def __init__(self, regions, default=0.0, geo_reference=None, verbose=False):
666
667        try:
668            len(regions)
669        except:
670            msg = 'Polygon_function takes a list of pairs (polygon, value).'
671            msg += 'Got %s' %regions
672            raise msg
673
674
675        T = regions[0]
676
677        if isinstance(T, basestring):
678            msg = 'You passed in a list of text values into polygon_function'
679            msg += ' instead of a list of pairs (polygon, value): "%s"' %T
680           
681            raise Exception, msg
682       
683        try:
684            a = len(T)
685        except:
686            msg = 'Polygon_function takes a list of pairs (polygon, value).'
687            msg += 'Got %s' %str(T)
688            raise msg
689
690        msg = 'Each entry in regions have two components: (polygon, value).'
691        msg +='I got %s' %str(T)
692        assert a == 2, msg
693
694
695        if geo_reference is None:
696            from anuga.coordinate_transforms.geo_reference import Geo_reference
697            geo_reference = Geo_reference()
698
699
700        self.default = default
701
702        # Make points in polygons relative to geo_reference
703        self.regions = []
704        for polygon, value in regions:
705            P = geo_reference.change_points_geo_ref(polygon)
706            self.regions.append( (P, value) )
707
708        self.verbose = verbose
709
710
711
712
713    def __call__(self, x, y):
714        x = numpy.array(x).astype(Float)
715        y = numpy.array(y).astype(Float)
716
717        assert len(x.shape) == 1 and len(y.shape) == 1
718       
719        N = x.shape[0]
720        assert y.shape[0] == N
721
722        points = numpy.ascontiguousarray(numpy.concatenate( (x[:,numpy.newaxis], y[:,numpy.newaxis]), axis=1 ))
723       
724        if callable(self.default):
725            z = self.default(x, y)
726        else:
727            z = numpy.ones(N, Float) * self.default
728
729        for polygon, value in self.regions:
730            indices = inside_polygon(points, polygon, verbose=self.verbose)
731
732            # FIXME: This needs to be vectorised
733            if callable(value):
734                for i in indices:
735                    xx = numpy.array([x[i]])
736                    yy = numpy.array([y[i]])
737                    z[i] = value(xx, yy)[0]
738            else:
739                for i in indices:
740                    z[i] = value
741
742        return z
743
744
745def read_polygon(filename, split=','):
746    """Read points assumed to form a polygon.
747       There must be exactly two numbers in each line separated by a comma.
748       No header.
749    """
750
751    #Get polygon
752    fid = open(filename)
753    lines = fid.readlines()
754    fid.close()
755    polygon = []
756    for line in lines:
757        fields = line.split(split)
758        polygon.append( [float(fields[0]), float(fields[1])] )
759
760    return polygon
761
762
763def write_polygon(polygon, filename=None):
764    """Write polygon to csv file.
765       There will be exactly two numbers, easting and northing,
766       in each line separated by a comma.
767       
768       No header.   
769    """
770
771    fid = open(filename, 'w')
772    for point in polygon:
773        fid.write('%f, %f\n' %point)
774    fid.close()
775   
776
777def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
778    """Populate given polygon with uniformly distributed points.
779
780    Input:
781       polygon - list of vertices of polygon
782       number_of_points - (optional) number of points
783       seed - seed for random number generator (default=None)
784       exclude - list of polygons (inside main polygon) from where points should be excluded
785
786    Output:
787       points - list of points inside polygon
788
789    Examples:
790       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
791       will return five randomly selected points inside the unit square
792    """
793
794    from random import uniform, seed as seed_function
795
796    seed_function(seed)
797
798    points = []
799
800    # Find outer extent of polygon
801    max_x = min_x = polygon[0][0]
802    max_y = min_y = polygon[0][1]
803    for point in polygon[1:]:
804        x = point[0]
805        if x > max_x: max_x = x
806        if x < min_x: min_x = x
807        y = point[1]
808        if y > max_y: max_y = y
809        if y < min_y: min_y = y
810
811
812    while len(points) < number_of_points:
813        x = uniform(min_x, max_x)
814        y = uniform(min_y, max_y)
815
816        append = False
817        if is_inside_polygon([x,y], polygon):
818
819            append = True
820
821            #Check exclusions
822            if exclude is not None:
823                for ex_poly in exclude:
824                    if is_inside_polygon([x,y], ex_poly):
825                        append = False
826
827
828        if append is True:
829            points.append([x,y])
830
831    return points
832
833
834def point_in_polygon(polygon, delta=1e-8):
835    """Return a point inside a given polygon which will be close to the
836    polygon edge.
837
838    Input:
839       polygon - list of vertices of polygon
840       delta - the square root of 2 * delta is the maximum distance from the
841       polygon points and the returned point.
842    Output:
843       points - a point inside polygon
844
845       searches in all diagonals and up and down (not left and right)
846    """
847    import exceptions
848    class Found(exceptions.Exception): pass
849
850    point_in = False
851    while not point_in:
852        try:
853            for poly_point in polygon: #[1:]:
854                for x_mult in range (-1,2):
855                    for y_mult in range (-1,2):
856                        x = poly_point[0]
857                        y = poly_point[1]
858                        if x == 0:
859                            x_delta = x_mult*delta
860                        else:
861                            x_delta = x+x_mult*x*delta
862
863                        if y == 0:
864                            y_delta = y_mult*delta
865                        else:
866                            y_delta = y+y_mult*y*delta
867
868                        point = [x_delta, y_delta]
869                        #print "point",point
870                        if is_inside_polygon(point, polygon, closed=False):
871                            raise Found
872        except Found:
873            point_in = True
874        else:
875            delta = delta*0.1
876    return point
877
878
879def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
880    """Calculate the approximate number of triangles inside the
881    bounding polygon and the other interior regions
882
883    Polygon areas are converted to square Kms
884
885    FIXME: Add tests for this function
886    """
887   
888    from anuga.utilities.polygon import polygon_area
889
890
891    # TO DO check if any of the regions fall inside one another
892
893    print '----------------------------------------------------------------------------'
894    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
895    print '----------------------------------------------------------------------------'   
896       
897    no_triangles = 0.0
898    area = polygon_area(bounding_poly)
899   
900    for poly, resolution in interior_regions:
901        this_area = polygon_area(poly)
902        this_triangles = this_area/resolution
903        no_triangles += this_triangles
904        area -= this_area
905       
906        print 'Interior ',
907        print ('%.0f' %resolution).ljust(25),
908        print ('%.2f' %(this_area/1000000)).ljust(19),
909        print '%d' %(this_triangles)
910       
911    bound_triangles = area/remainder_res
912    no_triangles += bound_triangles
913
914    print 'Bounding ',
915    print ('%.0f' %remainder_res).ljust(25),
916    print ('%.2f' %(area/1000000)).ljust(19),
917    print '%d' %(bound_triangles)   
918
919    total_number_of_triangles = no_triangles/0.7
920
921    print 'Estimated total number of triangles: %d' %total_number_of_triangles
922    print 'Note: This is generally about 20% less than the final amount'   
923
924    return int(total_number_of_triangles)
925
926
927def decimate_polygon(polygon, factor=10):
928    """Reduce number of points in polygon by the specified
929    factor (default=10, hence the name of the function) such that
930    the extrema in both axes are preserved.
931
932    Return reduced polygon
933    """
934
935    # FIXME(Ole): This doesn't work at present,
936    # but it isn't critical either
937
938    # Find outer extent of polygon
939    num_polygon = ensure_numeric(polygon)
940    max_x = max(num_polygon[:,0])
941    max_y = max(num_polygon[:,1])
942    min_x = min(num_polygon[:,0])
943    min_y = min(num_polygon[:,1])       
944
945    # Keep only some points making sure extrema are kept
946    reduced_polygon = []   
947    for i, point in enumerate(polygon):
948        x = point[0]
949        y = point[1]       
950        if x in [min_x, max_x] and y in [min_y, max_y]:
951            # Keep
952            reduced_polygon.append(point)
953        else:
954            if len(reduced_polygon)*factor < i:
955                reduced_polygon.append(point)               
956
957    return reduced_polygon
958
959##############################################
960#Initialise module
961
962from anuga.utilities.compile import can_use_C_extension
963if can_use_C_extension('polygon_ext.c'):
964    # Underlying C implementations can be accessed
965    from polygon_ext import _point_on_line
966    from polygon_ext import _separate_points_by_polygon
967    #from polygon_ext import _intersection
968
969else:
970    msg = 'C implementations could not be accessed by %s.\n ' %__file__
971    msg += 'Make sure compile_all.py has been run as described in '
972    msg += 'the ANUGA installation guide.'
973    raise Exception, msg
974
975
976if __name__ == "__main__":
977    pass
Note: See TracBrowser for help on using the repository browser.