# source:branches/numpy/anuga/utilities/polygon.py@6304

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

Initial commit of numpy changes. Still a long way to go.

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