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

Last change on this file since 8057 was 8057, checked in by steve, 9 years ago

Added a 'closed' argument to read_polygon so that the code correctly tests for a self intersecting polygon (both open and closed)

File size: 39.4 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[0], point[1],
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[0]/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[0]/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 line_intersect(triangles, line, verbose=False):
260    """Determine if a polygon and triangle overlap
261
262    """
263    line = ensure_numeric(line)
264    triangles = ensure_numeric(triangles)
265   
266    M = triangles.shape[0]/3  # Number of triangles
267
268    indices = num.zeros(M, num.int)
269
270    count = _line_intersect(line, triangles, indices)
271
272    if verbose:
273        log.critical('Found %d triangles (out of %d) that overlap the polygon' % (count, M))
274
275    return indices[:count]
276   
277def not_line_intersect(triangles, line, verbose=False):
278    """Determine if a polyline and triangle overlap
279
280    """
281    line = ensure_numeric(line)
282    triangles = ensure_numeric(triangles)
283   
284    M = triangles.shape[0]/3  # Number of triangles
285
286    indices = num.zeros(M, num.int)
287
288    count = _line_intersect(line, triangles, indices)
289
290    if verbose:
291        log.critical('Found %d triangles (out of %d) that intersect the line' % (count, M))
292
293    return indices[count:]   
294   
295
296def is_inside_triangle(point, triangle, 
297                       closed=True, 
298                       rtol=1.0e-12,
299                       atol=1.0e-12,                     
300                       check_inputs=True):
301    """Determine if one point is inside a triangle
302   
303    This uses the barycentric method:
304   
305    Triangle is A, B, C
306    Point P can then be written as
307   
308    P = A + alpha * (C-A) + beta * (B-A)
309    or if we let
310    v=P-A, v0=C-A, v1=B-A   
311   
312    v = alpha*v0 + beta*v1
313
314    Dot this equation by v0 and v1 to get two:
315   
316    dot(v0, v) = alpha*dot(v0, v0) + beta*dot(v0, v1)
317    dot(v1, v) = alpha*dot(v1, v0) + beta*dot(v1, v1)   
318   
319    or if a_ij = dot(v_i, v_j) and b_i = dot(v_i, v)
320    the matrix equation:
321   
322    a_00 a_01   alpha     b_0
323                       =
324    a_10 a_11   beta      b_1
325   
326    Solving for alpha and beta yields:
327   
328    alpha = (b_0*a_11 - b_1*a_01)/denom
329    beta =  (b_1*a_00 - b_0*a_10)/denom
330   
331    with denom = a_11*a_00 - a_10*a_01
332   
333    The point is in the triangle whenever
334    alpha and beta and their sums are in the unit interval.
335   
336    rtol and atol will determine how close the point has to be to the edge
337    before it is deemed to be on the edge.
338   
339    """
340
341    triangle = ensure_numeric(triangle)       
342    point = ensure_numeric(point, num.float)   
343   
344    if check_inputs is True:
345        msg = 'is_inside_triangle must be invoked with one point only'
346        assert num.allclose(point.shape, [2]), msg
347   
348   
349    # Use C-implementation
350    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
351   
352def is_complex(polygon, closed=True, verbose=False):
353    """Check if a polygon is complex (self-intersecting).
354       Uses a sweep algorithm that is O(n^2) in the worst case, but
355       for most normal looking polygons it'll be O(n log n).
356
357       polygon is a list of points that define a closed polygon.
358       verbose will print a list of the intersection points if true
359       
360       Return True if polygon is complex.
361    """           
362           
363    def key_xpos(item):
364        """ Return the x coord out of the passed point for sorting key. """
365        return (item[0][0])
366   
367    def segments_joined(seg0, seg1):
368        """ See if there are identical segments in the 2 lists. """
369        for i in seg0:
370            for j in seg1:   
371                if i == j: return True
372        return False
373       
374    polygon = ensure_numeric(polygon, num.float)
375
376    # build a list of discrete segments from the polygon
377    unsorted_segs = []
378    for i in range(0, len(polygon)-1):
379        unsorted_segs.append([list(polygon[i]), list(polygon[i+1])])
380
381    if closed:
382        unsorted_segs.append([list(polygon[0]), list(polygon[-1])])
383   
384    # all segments must point in same direction
385    for val in unsorted_segs:
386        if val[0][0] > val[1][0]:
387            val[0], val[1] = val[1], val[0]   
388           
389    l_x = sorted(unsorted_segs, key=key_xpos)
390
391    comparisons = 0
392   
393    # loop through, only comparing lines that partially overlap in x
394    for index, leftmost in enumerate(l_x):
395        cmp = index+1
396        while cmp < len(l_x) and leftmost[1][0] > l_x[cmp][0][0]:
397            if not segments_joined(leftmost, l_x[cmp]):
398                (type, point) = intersection(leftmost, l_x[cmp])
399                comparisons += 1
400                if type != 0 and type != 4 or (type == 2 and list(point[0]) !=\
401                                                list(point[1])):
402                    if verbose:
403                        print 'Self-intersecting polygon found, type ', type
404                        print 'point', point,
405                        print 'vertices: ', leftmost, ' - ', l_x[cmp] 
406                    return True           
407            cmp += 1
408       
409    return False
410   
411
412def is_inside_polygon(point, polygon, closed=True, verbose=False):
413    """Determine if one point is inside a polygon
414
415    See inside_polygon for more details
416    """
417
418    indices = inside_polygon(point, polygon, closed, verbose)
419
420    if indices.shape[0] == 1:
421        return True
422    elif indices.shape[0] == 0:
423        return False
424    else:
425        msg = 'is_inside_polygon must be invoked with one point only'
426        raise Exception(msg)
427
428##
429# @brief Determine which of a set of points are inside a polygon.
430# @param points A set of points (tuple, list or array).
431# @param polygon A set of points defining a polygon (tuple, list or array).
432# @param closed True if points on boundary are considered 'inside' polygon.
433# @param verbose True if this function is to be verbose.
434# @return A list of indices of points inside the polygon.
435def inside_polygon(points, polygon, closed=True, verbose=False):
436    """Determine points inside a polygon
437
438       Functions inside_polygon and outside_polygon have been defined in
439       terms of separate_by_polygon which will put all inside indices in
440       the first part of the indices array and outside indices in the last
441
442       See separate_points_by_polygon for documentation
443
444       points and polygon can be a geospatial instance,
445       a list or a numeric array
446    """
447
448    try:
449        points = ensure_absolute(points)
450    except NameError, err:
451        raise NameError, err
452    except:
453        # If this fails it is going to be because the points can't be
454        # converted to a numeric array.
455        msg = 'Points could not be converted to numeric array' 
456        raise Exception, msg
457
458    polygon = ensure_absolute(polygon)       
459    try:
460        polygon = ensure_absolute(polygon)
461    except NameError, e:
462        raise NameError, e
463    except:
464        # If this fails it is going to be because the points can't be
465        # converted to a numeric array.
466        msg = ('Polygon %s could not be converted to numeric array'
467               % (str(polygon)))
468        raise Exception, msg
469
470    if len(points.shape) == 1:
471        # Only one point was passed in. Convert to array of points
472        points = num.reshape(points, (1,2))
473
474    indices, count = separate_points_by_polygon(points, polygon,
475                                                closed=closed,
476                                                verbose=verbose)
477
478    # Return indices of points inside polygon
479    return indices[:count]
480
481##
482# @brief Determine if one point is outside a polygon.
483# @param point The point of interest.
484# @param polygon The polygon to test inclusion in.
485# @param closed True if points on boundary are considered 'inside' polygon.
486# @param verbose True if this function is to be verbose.
487# @return True if point is outside the polygon.
488# @note Uses inside_polygon() to do the work.
489def is_outside_polygon(point, polygon, closed=True, verbose=False,
490                       points_geo_ref=None, polygon_geo_ref=None):
491    """Determine if one point is outside a polygon
492
493    See outside_polygon for more details
494    """
495
496    indices = outside_polygon(point, polygon, closed, verbose)
497
498    if indices.shape[0] == 1:
499        return True
500    elif indices.shape[0] == 0:
501        return False
502    else:
503        msg = 'is_outside_polygon must be invoked with one point only'
504        raise Exception, msg
505
506##
507# @brief Determine which of a set of points are outside a polygon.
508# @param points A set of points (tuple, list or array).
509# @param polygon A set of points defining a polygon (tuple, list or array).
510# @param closed True if points on boundary are considered 'inside' polygon.
511# @param verbose True if this function is to be verbose.
512# @return A list of indices of points outside the polygon.
513def outside_polygon(points, polygon, closed = True, verbose = False):
514    """Determine points outside a polygon
515
516       Functions inside_polygon and outside_polygon have been defined in
517       terms of separate_by_polygon which will put all inside indices in
518       the first part of the indices array and outside indices in the last
519
520       See separate_points_by_polygon for documentation
521    """
522
523    try:
524        points = ensure_numeric(points, num.float)
525    except NameError, e:
526        raise NameError, e
527    except:
528        msg = 'Points could not be converted to numeric array'
529        raise Exception, msg
530
531    try:
532        polygon = ensure_numeric(polygon, num.float)
533    except NameError, e:
534        raise NameError, e
535    except:
536        msg = 'Polygon could not be converted to numeric array'
537        raise Exception, msg
538
539    if len(points.shape) == 1:
540        # Only one point was passed in. Convert to array of points
541        points = num.reshape(points, (1, 2))
542
543    indices, count = separate_points_by_polygon(points, polygon,
544                                                closed=closed,
545                                                verbose=verbose)
546
547    # Return indices of points outside polygon
548    if count == len(indices):
549        # No points are outside
550        return num.array([])
551    else:
552        return indices[count:][::-1]  #return reversed
553
554##
555# @brief Separate a list of points into two sets inside+outside a polygon.
556# @param points A set of points (tuple, list or array).
557# @param polygon A set of points defining a polygon (tuple, list or array).
558# @param closed True if points on boundary are considered 'inside' polygon.
559# @param verbose True if this function is to be verbose.
560# @return A tuple (in, out) of point indices for poinst inside amd outside.
561def in_and_outside_polygon(points, polygon, closed=True, verbose=False):
562    """Determine points inside and outside a polygon
563
564       See separate_points_by_polygon for documentation
565
566       Returns an array of points inside and array of points outside the polygon
567    """
568
569    try:
570        points = ensure_numeric(points, num.float)
571    except NameError, e:
572        raise NameError, e
573    except:
574        msg = 'Points could not be converted to numeric array'
575        raise Exception, msg
576
577    try:
578        polygon = ensure_numeric(polygon, num.float)
579    except NameError, e:
580        raise NameError, e
581    except:
582        msg = 'Polygon could not be converted to numeric array'
583        raise Exception, msg
584
585    if len(points.shape) == 1:
586        # Only one point was passed in. Convert to array of points
587        points = num.reshape(points, (1, 2))
588
589    indices, count = separate_points_by_polygon(points, polygon,
590                                                closed=closed,
591                                                verbose=verbose)
592
593    # Returns indices of points inside and indices of points outside
594    # the polygon
595    if count == len(indices):
596        # No points are outside
597        return indices[:count], []
598    else:
599        return  indices[:count], indices[count:][::-1]  #return reversed
600
601
602
603def separate_points_by_polygon(points, polygon,
604                               closed=True, 
605                               check_input=True,
606                               verbose=False):
607    """Determine whether points are inside or outside a polygon
608
609    Input:
610       points - Tuple of (x, y) coordinates, or list of tuples
611       polygon - list of vertices of polygon
612       closed - (optional) determine whether points on boundary should be
613       regarded as belonging to the polygon (closed = True)
614       or not (closed = False)
615       check_input: Allows faster execution if set to False
616
617    Outputs:
618       indices: array of same length as points with indices of points falling
619       inside the polygon listed from the beginning and indices of points
620       falling outside listed from the end.
621
622       count: count of points falling inside the polygon
623
624       The indices of points inside are obtained as indices[:count]
625       The indices of points outside are obtained as indices[count:]
626
627    Examples:
628       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
629
630       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
631       will return the indices [0, 2, 1] and count == 2 as only the first
632       and the last point are inside the unit square
633
634    Remarks:
635       The vertices may be listed clockwise or counterclockwise and
636       the first point may optionally be repeated.
637       Polygons do not need to be convex.
638       Polygons can have holes in them and points inside a hole is
639       regarded as being outside the polygon.
640
641    Algorithm is based on work by Darel Finley,
642    http://www.alienryderflex.com/polygon/
643
644    Uses underlying C-implementation in polygon_ext.c
645    """
646
647    if check_input:
648        #Input checks
649        assert isinstance(closed, bool), \
650                    'Keyword argument "closed" must be boolean'
651        assert isinstance(verbose, bool), \
652                    'Keyword argument "verbose" must be boolean'
653
654        try:
655            points = ensure_numeric(points, num.float)
656        except NameError, e:
657            raise NameError, e
658        except:
659            msg = 'Points could not be converted to numeric array'
660            raise Exception(msg)
661
662        try:
663            polygon = ensure_numeric(polygon, num.float)
664        except NameError, e:
665            raise NameError(e)
666        except:
667            msg = 'Polygon could not be converted to numeric array'
668            raise Exception(msg)
669
670        msg = 'Polygon array must be a 2d array of vertices'
671        assert len(polygon.shape) == 2, msg
672
673        msg = 'Polygon array must have two columns' 
674        assert polygon.shape[1] == 2, msg
675
676
677        msg = ('Points array must be 1 or 2 dimensional. '
678               'I got %d dimensions' % len(points.shape))
679        assert 0 < len(points.shape) < 3, msg
680
681        if len(points.shape) == 1:
682            # Only one point was passed in.  Convert to array of points.
683            points = num.reshape(points, (1, 2))
684   
685            msg = ('Point array must have two columns (x,y), '
686                   'I got points.shape[1]=%d' % points.shape[0])
687            assert points.shape[1]==2, msg
688
689       
690            msg = ('Points array must be a 2d array. I got %s.'
691                   % str(points[:30]))
692            assert len(points.shape) == 2, msg
693
694            msg = 'Points array must have two columns'
695            assert points.shape[1] == 2, msg
696
697    N = polygon.shape[0] # Number of vertices in polygon
698    M = points.shape[0]  # Number of points
699
700    indices = num.zeros(M, num.int)
701
702    count = _separate_points_by_polygon(points, polygon, indices,
703                                        int(closed), int(verbose))
704
705    if verbose:
706        log.critical('Found %d points (out of %d) inside polygon' % (count, M))
707
708    return indices, count
709
710
711def polygon_area(input_polygon):
712    """ Determine area of arbitrary polygon.
713
714        input_polygon The polygon to get area of.
715       
716        return A scalar value for the polygon area.
717
718        Reference:     http://mathworld.wolfram.com/PolygonArea.html
719    """
720    # Move polygon to origin (0,0) to avoid rounding errors
721    # This makes a copy of the polygon to avoid destroying it
722    input_polygon = ensure_numeric(input_polygon)
723    min_x = min(input_polygon[:, 0])
724    min_y = min(input_polygon[:, 1])
725    polygon = input_polygon - [min_x, min_y]
726
727    # Compute area
728    n = len(polygon)
729    poly_area = 0.0
730
731    for i in range(n):
732        pti = polygon[i]
733        if i == n-1:
734            pt1 = polygon[0]
735        else:
736            pt1 = polygon[i+1]
737        xi = pti[0]
738        yi1 = pt1[1]
739        xi1 = pt1[0]
740        yi = pti[1]
741        poly_area += xi*yi1 - xi1*yi
742
743    return abs(poly_area/2)
744
745
746def plot_polygons(polygons_points,
747                  style=None,
748                  figname=None,
749                  label=None,
750                  alpha=None):
751    """ Take list of polygons and plot.
752
753    Inputs:
754
755    polygons         - list of polygons
756
757    style            - style list corresponding to each polygon
758                     - for a polygon, use 'line'
759                     - for points falling outside a polygon, use 'outside'
760                     - style can also be user defined as in normal pylab plot.
761
762    figname          - name to save figure to
763
764    label            - title for plotA
765
766    alpha            - transparency of polygon fill, 0.0=none, 1.0=solid
767                       if not supplied, no fill.
768
769    Outputs:
770
771    - plot of polygons
772    """
773
774    from pylab import ion, hold, plot, savefig, xlabel, \
775                      ylabel, title, close, title, fill
776
777    assert type(polygons_points) == list, \
778                'input must be a list of polygons and/or points'
779
780    ion()
781    hold(True)
782
783    if label is None:
784        label = ''
785
786    # clamp alpha to sensible range
787    if alpha:
788        try:
789            alpha = float(alpha)
790        except ValueError:
791            alpha = None
792        else:
793            alpha = max(0.0, min(1.0, alpha))
794
795    num_points = len(polygons_points)
796    colour = []
797    if style is None:
798        style_type = 'line'
799        style = []
800        for i in range(num_points):
801            style.append(style_type)
802            colour.append('b-')
803    else:
804        for style_name in style:
805            if style_name == 'line':
806                colour.append('b-')
807            if style_name == 'outside':
808                colour.append('r.')
809            if style_name == 'point':
810                colour.append('g.')
811            if style_name not in ['line', 'outside', 'point']:
812                colour.append(style_name)
813
814    for i, item in enumerate(polygons_points):
815        pt_x, pt_y = _poly_xy(item)
816        plot(pt_x, pt_y, colour[i])
817        if alpha:
818            fill(pt_x, pt_y, colour[i], alpha=alpha)
819        xlabel('x')
820        ylabel('y')
821        title(label)
822
823    if figname is not None:
824        savefig(figname)
825    else:
826        savefig('test_image')
827
828    close('all')
829
830
831def _poly_xy(polygon):
832    """ this is used within plot_polygons so need to duplicate
833        the first point so can have closed polygon in plot
834        # @param polygon A set of points defining a polygon.
835        # @param verbose True if this function is to be verbose.
836        # @return A tuple (x, y) of X and Y coordinates of the polygon.
837        # @note We duplicate the first point so can have closed polygon in plot.
838    """
839
840    try:
841        polygon = ensure_numeric(polygon, num.float)
842    except NameError, err:
843        raise NameError, err
844    except:
845        msg = ('Polygon %s could not be converted to numeric array'
846               % (str(polygon)))
847        raise Exception, msg
848
849    pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0)
850    pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0)
851
852    return pts_x, pts_y
853
854
855################################################################################
856# Functions to read and write polygon information
857################################################################################
858
859def read_polygon(filename, delimiter=',', closed='True'):
860    """ Read points assumed to form a (closed) polygon.
861        Can also be used toread  in a polyline (closed = 'False')
862
863        Also checks to make sure polygon (polyline)
864        is not complex (self-intersecting).
865
866        filename Path to file containing polygon data.
867        delimiter Delimiter to split polygon data with.
868        A list of point data from the polygon file.
869
870        There must be exactly two numbers in each line separated by the delimiter.
871        No header.
872    """
873
874    fid = open(filename)
875    lines = fid.readlines()
876    fid.close()
877    polygon = []
878    for line in lines:
879        fields = line.split(delimiter)
880        polygon.append([float(fields[0]), float(fields[1])])
881   
882    # check this is a valid polygon (polyline).
883    if is_complex(polygon, closed, verbose=True):
884        msg = 'ERROR: Self-intersecting polygon detected in file '
885        msg += filename +'. A complex polygon will not '
886        msg += 'necessarily break the algorithms within ANUGA, but it'
887        msg += 'usually signifies pathological data. Please fix this file.'
888        raise Exception, msg
889   
890    return polygon
891
892
893def write_polygon(polygon, filename=None):
894    """Write polygon to csv file.
895
896    There will be exactly two numbers, easting and northing, in each line
897    separated by a comma.
898
899    No header.
900    """
901
902    fid = open(filename, 'w')
903    for point in polygon:
904        fid.write('%f, %f\n' % point)
905    fid.close()
906
907
908def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
909    """Populate given polygon with uniformly distributed points.
910
911    Input:
912       polygon - list of vertices of polygon
913       number_of_points - (optional) number of points
914       seed - seed for random number generator (default=None)
915       exclude - list of polygons (inside main polygon) from where points
916                 should be excluded
917
918    Output:
919       points - list of points inside polygon
920
921    Examples:
922       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
923       will return five randomly selected points inside the unit square
924    """
925
926    from random import uniform, seed as seed_function
927
928    seed_function(seed)
929
930    points = []
931
932    # Find outer extent of polygon
933    extents = AABB(polygon)
934   
935    while len(points) < number_of_points:
936        rand_x = uniform(extents.xmin, extents.xmax)
937        rand_y = uniform(extents.ymin, extents.ymax)
938
939        append = False
940        if is_inside_polygon([rand_x, rand_y], polygon):
941            append = True
942
943            #Check exclusions
944            if exclude is not None:
945                for ex_poly in exclude:
946                    if is_inside_polygon([rand_x, rand_y], ex_poly):
947                        append = False
948
949        if append is True:
950            points.append([rand_x, rand_y])
951
952    return points
953
954
955def point_in_polygon(polygon, delta=1e-8):
956    """Return a point inside a given polygon which will be close to the
957    polygon edge.
958
959    Input:
960       polygon - list of vertices of polygon
961       delta - the square root of 2 * delta is the maximum distance from the
962       polygon points and the returned point.
963    Output:
964       points - a point inside polygon
965
966       searches in all diagonals and up and down (not left and right).
967    """
968
969    polygon = ensure_numeric(polygon)
970   
971    while True:
972        for poly_point in polygon:
973            for x_mult in range(-1, 2):
974                for y_mult in range(-1, 2):
975                    pt_x, pt_y = poly_point
976
977                    if pt_x == 0:
978                        x_delta = x_mult * delta
979                    else:
980                        x_delta = pt_x + x_mult*pt_x*delta
981
982                    if pt_y == 0:
983                        y_delta = y_mult * delta
984                    else:
985                        y_delta = pt_y + y_mult*pt_y*delta
986
987                    point = [x_delta, y_delta]
988
989                    if is_inside_polygon(point, polygon, closed=False):
990                        return point
991        delta = delta * 0.1
992
993
994def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
995    """Calculate the approximate number of triangles inside the
996    bounding polygon and the other interior regions
997
998    Polygon areas are converted to square Kms
999
1000    FIXME: Add tests for this function
1001    """
1002
1003    # TO DO check if any of the regions fall inside one another
1004
1005    log.critical('-' * 80)
1006    log.critical('Polygon   Max triangle area (m^2)   Total area (km^2) '
1007                 'Estimated #triangles')
1008    log.critical('-' * 80)
1009       
1010    no_triangles = 0.0
1011    area = polygon_area(bounding_poly)
1012
1013    for poly, resolution in interior_regions:
1014        this_area = polygon_area(poly)
1015        this_triangles = this_area/resolution
1016        no_triangles += this_triangles
1017        area -= this_area
1018
1019        log.critical('Interior %s%s%d'
1020                     % (('%.0f' % resolution).ljust(25),
1021                        ('%.2f' % (this_area/1000000)).ljust(19), 
1022                        this_triangles))
1023        #print 'Interior ',
1024        #print ('%.0f' % resolution).ljust(25),
1025        #print ('%.2f' % (this_area/1000000)).ljust(19),
1026        #print '%d' % (this_triangles)
1027
1028    bound_triangles = area/remainder_res
1029    no_triangles += bound_triangles
1030
1031    log.critical('Bounding %s%s%d'
1032                 % (('%.0f' % remainder_res).ljust(25),
1033                    ('%.2f' % (area/1000000)).ljust(19),
1034                    bound_triangles))
1035    #print 'Bounding ',
1036    #print ('%.0f' % remainder_res).ljust(25),
1037    #print ('%.2f' % (area/1000000)).ljust(19),
1038    #print '%d' % (bound_triangles)
1039
1040    total_number_of_triangles = no_triangles/0.7
1041
1042    log.critical('Estimated total number of triangles: %d'
1043                 % total_number_of_triangles)
1044    log.critical('Note: This is generally about 20%% '
1045                 'less than the final amount')
1046
1047    return int(total_number_of_triangles)
1048
1049
1050def decimate_polygon(polygon, factor=10):
1051    """Reduce number of points in polygon by the specified
1052    factor (default=10, hence the name of the function) such that
1053    the extrema in both axes are preserved.
1054
1055##
1056# @brief Reduce number of points in polygon by the specified factor.
1057# @param polygon The polygon to reduce.
1058# @param factor The factor to reduce polygon points by (default 10).
1059# @note The extrema of both axes are preserved.
1060
1061    Return reduced polygon
1062    """
1063
1064    # FIXME(Ole): This doesn't work at present,
1065    # but it isn't critical either
1066
1067    # Find outer extent of polygon
1068    num_polygon = ensure_numeric(polygon)
1069    max_x = max(num_polygon[:, 0])
1070    max_y = max(num_polygon[:, 1])
1071    min_x = min(num_polygon[:, 0])
1072    min_y = min(num_polygon[:, 1])
1073
1074    # Keep only some points making sure extrema are kept
1075    reduced_polygon = []
1076    for i, point in enumerate(polygon):
1077        if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]:
1078            # Keep
1079            reduced_polygon.append(point)
1080        else:
1081            if len(reduced_polygon)*factor < i:
1082                reduced_polygon.append(point)
1083
1084    return reduced_polygon
1085
1086
1087def interpolate_polyline(data,
1088                         polyline_nodes,
1089                         gauge_neighbour_id,
1090                         interpolation_points=None,
1091                         rtol=1.0e-6,
1092                         atol=1.0e-8):
1093    """Interpolate linearly between values data on polyline nodes
1094    of a polyline to list of interpolation points.
1095
1096    data is the data on the polyline nodes.
1097
1098    Inputs:
1099      data: Vector or array of data at the polyline nodes.
1100      polyline_nodes: Location of nodes where data is available.
1101      gauge_neighbour_id: ?
1102      interpolation_points: Interpolate polyline data to these positions.
1103          List of coordinate pairs [x, y] of
1104          data points or an nx2 numeric array or a Geospatial_data object
1105      rtol, atol: Used to determine whether a point is on the polyline or not.
1106                  See point_on_line.
1107
1108    Output:
1109      Interpolated values at interpolation points
1110    """
1111
1112    if isinstance(interpolation_points, Geospatial_data):
1113        interpolation_points = interpolation_points.\
1114                                    get_data_points(absolute=True)
1115
1116    interpolated_values = num.zeros(len(interpolation_points), num.float)
1117
1118    data = ensure_numeric(data, num.float)
1119    polyline_nodes = ensure_numeric(polyline_nodes, num.float)
1120    interpolation_points = ensure_numeric(interpolation_points, num.float)
1121    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
1122
1123    num_nodes = polyline_nodes.shape[0]    # Number of nodes in polyline
1124
1125    # Input sanity check
1126    assert_msg = 'interpolation_points are not given (interpolate.py)'
1127    assert interpolation_points is not None, assert_msg
1128
1129    assert_msg = 'function value must be specified at every interpolation node'
1130    assert data.shape[0] == polyline_nodes.shape[0], assert_msg
1131
1132    assert_msg = 'Must define function value at one or more nodes'
1133    assert data.shape[0] > 0, assert_msg
1134
1135    if num_nodes == 1:
1136        assert_msg = 'Polyline contained only one point. I need more. '
1137        assert_msg += str(data)
1138        raise Exception, assert_msg
1139    elif num_nodes > 1:
1140        _interpolate_polyline(data,
1141                              polyline_nodes,
1142                              gauge_neighbour_id,
1143                              interpolation_points,
1144                              interpolated_values,
1145                              rtol,
1146                              atol)
1147
1148
1149    return interpolated_values
1150
1151   
1152def polylist2points_verts(polylist):
1153    """ Convert a list of polygons to discrete points and vertices.
1154    """
1155   
1156    offset = 0
1157    points = []
1158    vertices = []
1159    for poly in polylist:
1160        points.extend(poly)
1161        vertices.extend([[i, i+1] for i in range(offset, offset+len(poly)-1)])
1162        offset += len(poly)
1163               
1164    return points, vertices
1165
1166
1167################################################################################
1168# Initialise module
1169################################################################################
1170
1171from anuga.utilities import compile
1172if compile.can_use_C_extension('polygon_ext.c'):
1173    # Underlying C implementations can be accessed
1174    from polygon_ext import _point_on_line
1175    from polygon_ext import _separate_points_by_polygon
1176    from polygon_ext import _interpolate_polyline   
1177    from polygon_ext import _polygon_overlap
1178    from polygon_ext import _line_intersect
1179    from polygon_ext import _is_inside_triangle       
1180    #from polygon_ext import _intersection
1181
1182else:
1183    ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__
1184    ERROR_MSG += 'Make sure compile_all.py has been run as described in '
1185    ERROR_MSG += 'the ANUGA installation guide.'
1186    raise Exception(ERROR_MSG)
1187
1188
1189if __name__ == "__main__":
1190    pass
Note: See TracBrowser for help on using the repository browser.