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

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

Cleaned up unit tests to use API.

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