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

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

Refactorings to increase code quality, fixed missing log import from sww_interrogate.

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