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

Last change on this file since 6223 was 6223, checked in by ole, 14 years ago

Various cleanup

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