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

Last change on this file since 7841 was 7841, checked in by hudson, 13 years ago

Refactorings to allow tests to pass.

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