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

Last change on this file since 6116 was 6116, checked in by ole, 15 years ago

Feeble start on reader of multiple polygons using csv2dict

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