# source:anuga_core/source/anuga/utilities/polygon.py@7686

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

Tickets 346/257: Added a simple brute-force is_complex function to find complex (ie pathological/self intersecting) polygons.
Added unit test to test above.

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