# source:trunk/anuga_core/source/anuga/geometry/polygon.py@8031

Last change on this file since 8031 was 8031, checked in by habili, 13 years ago

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