source: trunk/anuga_core/source/anuga/geometry/polygon.py @ 8151

Last change on this file since 8151 was 8151, checked in by wilsonr, 14 years ago

Removed '@brief' comments.

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