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

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

NumPy? conversion.

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