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

Last change on this file since 8021 was 8009, checked in by steve, 13 years ago

Slight conflict with init.py due to addition of hole_tags

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, 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):
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, ), 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)
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), list(polygon[-1])])
308
309    # all segments must point in same direction
310    for val in unsorted_segs:
311        if val > val:
312            val, val = val, val
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 > l_x[cmp]:
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) !=\
326                                                list(point)):
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 == 1:
346        return True
347    elif indices.shape == 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 == 1:
424        return True
425    elif indices.shape == 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 == 2, msg
600
601
602        msg = ('Points array must be 1 or 2 dimensional. '
603               'I got %d dimensions' % len(points.shape))
604        assert 0 < len(points.shape) < 3, msg
605
606        if len(points.shape) == 1:
607            # Only one point was passed in.  Convert to array of points.
608            points = num.reshape(points, (1, 2))
609
610            msg = ('Point array must have two columns (x,y), '
611                   'I got points.shape=%d' % points.shape)
612            assert points.shape==2, msg
613
614
615            msg = ('Points array must be a 2d array. I got %s.'
616                   % str(points[:30]))
617            assert len(points.shape) == 2, msg
618
619            msg = 'Points array must have two columns'
620            assert points.shape == 2, msg
621
622    N = polygon.shape # Number of vertices in polygon
623    M = points.shape  # Number of points
624
625    indices = num.zeros(M, num.int)
626
627    count = _separate_points_by_polygon(points, polygon, indices,
628                                        int(closed), int(verbose))
629
630    if verbose:
631        log.critical('Found %d points (out of %d) inside polygon' % (count, M))
632
633    return indices, count
634
635
636def polygon_area(input_polygon):
637    """ Determine area of arbitrary polygon.
638
639        input_polygon The polygon to get area of.
640
641        return A scalar value for the polygon area.
642
643        Reference:     http://mathworld.wolfram.com/PolygonArea.html
644    """
645    # Move polygon to origin (0,0) to avoid rounding errors
646    # This makes a copy of the polygon to avoid destroying it
647    input_polygon = ensure_numeric(input_polygon)
648    min_x = min(input_polygon[:, 0])
649    min_y = min(input_polygon[:, 1])
650    polygon = input_polygon - [min_x, min_y]
651
652    # Compute area
653    n = len(polygon)
654    poly_area = 0.0
655
656    for i in range(n):
657        pti = polygon[i]
658        if i == n-1:
659            pt1 = polygon
660        else:
661            pt1 = polygon[i+1]
662        xi = pti
663        yi1 = pt1
664        xi1 = pt1
665        yi = pti
666        poly_area += xi*yi1 - xi1*yi
667
668    return abs(poly_area/2)
669
670
671def plot_polygons(polygons_points,
672                  style=None,
673                  figname=None,
674                  label=None,
675                  alpha=None):
676    """ Take list of polygons and plot.
677
678    Inputs:
679
680    polygons         - list of polygons
681
682    style            - style list corresponding to each polygon
683                     - for a polygon, use 'line'
684                     - for points falling outside a polygon, use 'outside'
685                     - style can also be user defined as in normal pylab plot.
686
687    figname          - name to save figure to
688
689    label            - title for plotA
690
691    alpha            - transparency of polygon fill, 0.0=none, 1.0=solid
692                       if not supplied, no fill.
693
694    Outputs:
695
696    - plot of polygons
697    """
698
699    from pylab import ion, hold, plot, savefig, xlabel, \
700                      ylabel, title, close, title, fill
701
702    assert type(polygons_points) == list, \
703                'input must be a list of polygons and/or points'
704
705    ion()
706    hold(True)
707
708    if label is None:
709        label = ''
710
711    # clamp alpha to sensible range
712    if alpha:
713        try:
714            alpha = float(alpha)
715        except ValueError:
716            alpha = None
717        else:
718            alpha = max(0.0, min(1.0, alpha))
719
720    num_points = len(polygons_points)
721    colour = []
722    if style is None:
723        style_type = 'line'
724        style = []
725        for i in range(num_points):
726            style.append(style_type)
727            colour.append('b-')
728    else:
729        for style_name in style:
730            if style_name == 'line':
731                colour.append('b-')
732            if style_name == 'outside':
733                colour.append('r.')
734            if style_name == 'point':
735                colour.append('g.')
736            if style_name not in ['line', 'outside', 'point']:
737                colour.append(style_name)
738
739    for i, item in enumerate(polygons_points):
740        pt_x, pt_y = _poly_xy(item)
741        plot(pt_x, pt_y, colour[i])
742        if alpha:
743            fill(pt_x, pt_y, colour[i], alpha=alpha)
744        xlabel('x')
745        ylabel('y')
746        title(label)
747
748    if figname is not None:
749        savefig(figname)
750    else:
751        savefig('test_image')
752
753    close('all')
754
755
756def _poly_xy(polygon):
757    """ this is used within plot_polygons so need to duplicate
758        the first point so can have closed polygon in plot
759        # @param polygon A set of points defining a polygon.
760        # @param verbose True if this function is to be verbose.
761        # @return A tuple (x, y) of X and Y coordinates of the polygon.
762        # @note We duplicate the first point so can have closed polygon in plot.
763    """
764
765    try:
766        polygon = ensure_numeric(polygon, num.float)
767    except NameError, err:
768        raise NameError, err
769    except:
770        msg = ('Polygon %s could not be converted to numeric array'
771               % (str(polygon)))
772        raise Exception, msg
773
774    pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0)
775    pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0)
776
777    return pts_x, pts_y
778
779
780################################################################################
781# Functions to read and write polygon information
782################################################################################
783
785    """ Read points assumed to form a polygon.
786
787        Also checks to make sure polygon is not complex (self-intersecting).
788
789        filename Path to file containing polygon data.
790        delimiter Delimiter to split polygon data with.
791        A list of point data from the polygon file.
792
793        There must be exactly two numbers in each line separated by the delimiter.
795    """
796
797    fid = open(filename)
799    fid.close()
800    polygon = []
801    for line in lines:
802        fields = line.split(delimiter)
803        polygon.append([float(fields), float(fields)])
804
805    # check this is a valid polygon.
806    if is_complex(polygon, verbose=True):
807        msg = 'ERROR: Self-intersecting polygon detected in file '
808        msg += filename +'. A complex polygon will not '
809        msg += 'necessarily break the algorithms within ANUGA, but it'
810        msg += 'usually signifies pathological data. Please fix this file.'
811        raise Exception, msg
812
813    return polygon
814
815
816def write_polygon(polygon, filename=None):
817    """Write polygon to csv file.
818
819    There will be exactly two numbers, easting and northing, in each line
820    separated by a comma.
821
823    """
824
825    fid = open(filename, 'w')
826    for point in polygon:
827        fid.write('%f, %f\n' % point)
828    fid.close()
829
830
831def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
832    """Populate given polygon with uniformly distributed points.
833
834    Input:
835       polygon - list of vertices of polygon
836       number_of_points - (optional) number of points
837       seed - seed for random number generator (default=None)
838       exclude - list of polygons (inside main polygon) from where points
839                 should be excluded
840
841    Output:
842       points - list of points inside polygon
843
844    Examples:
845       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
846       will return five randomly selected points inside the unit square
847    """
848
849    from random import uniform, seed as seed_function
850
851    seed_function(seed)
852
853    points = []
854
855    # Find outer extent of polygon
856    extents = AABB(polygon)
857
858    while len(points) < number_of_points:
859        rand_x = uniform(extents.xmin, extents.xmax)
860        rand_y = uniform(extents.ymin, extents.ymax)
861
862        append = False
863        if is_inside_polygon([rand_x, rand_y], polygon):
864            append = True
865
866            #Check exclusions
867            if exclude is not None:
868                for ex_poly in exclude:
869                    if is_inside_polygon([rand_x, rand_y], ex_poly):
870                        append = False
871
872        if append is True:
873            points.append([rand_x, rand_y])
874
875    return points
876
877
878def point_in_polygon(polygon, delta=1e-8):
879    """Return a point inside a given polygon which will be close to the
880    polygon edge.
881
882    Input:
883       polygon - list of vertices of polygon
884       delta - the square root of 2 * delta is the maximum distance from the
885       polygon points and the returned point.
886    Output:
887       points - a point inside polygon
888
889       searches in all diagonals and up and down (not left and right).
890    """
891
892    polygon = ensure_numeric(polygon)
893
894    while True:
895        for poly_point in polygon:
896            for x_mult in range(-1, 2):
897                for y_mult in range(-1, 2):
898                    pt_x, pt_y = poly_point
899
900                    if pt_x == 0:
901                        x_delta = x_mult * delta
902                    else:
903                        x_delta = pt_x + x_mult*pt_x*delta
904
905                    if pt_y == 0:
906                        y_delta = y_mult * delta
907                    else:
908                        y_delta = pt_y + y_mult*pt_y*delta
909
910                    point = [x_delta, y_delta]
911
912                    if is_inside_polygon(point, polygon, closed=False):
913                        return point
914        delta = delta * 0.1
915
916
917def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
918    """Calculate the approximate number of triangles inside the
919    bounding polygon and the other interior regions
920
921    Polygon areas are converted to square Kms
922
923    FIXME: Add tests for this function
924    """
925
926    # TO DO check if any of the regions fall inside one another
927
928    log.critical('-' * 80)
929    log.critical('Polygon   Max triangle area (m^2)   Total area (km^2) '
930                 'Estimated #triangles')
931    log.critical('-' * 80)
932
933    no_triangles = 0.0
934    area = polygon_area(bounding_poly)
935
936    for poly, resolution in interior_regions:
937        this_area = polygon_area(poly)
938        this_triangles = this_area/resolution
939        no_triangles += this_triangles
940        area -= this_area
941
942        log.critical('Interior %s%s%d'
943                     % (('%.0f' % resolution).ljust(25),
944                        ('%.2f' % (this_area/1000000)).ljust(19),
945                        this_triangles))
946        #print 'Interior ',
947        #print ('%.0f' % resolution).ljust(25),
948        #print ('%.2f' % (this_area/1000000)).ljust(19),
949        #print '%d' % (this_triangles)
950
951    bound_triangles = area/remainder_res
952    no_triangles += bound_triangles
953
954    log.critical('Bounding %s%s%d'
955                 % (('%.0f' % remainder_res).ljust(25),
956                    ('%.2f' % (area/1000000)).ljust(19),
957                    bound_triangles))
958    #print 'Bounding ',
959    #print ('%.0f' % remainder_res).ljust(25),
960    #print ('%.2f' % (area/1000000)).ljust(19),
961    #print '%d' % (bound_triangles)
962
963    total_number_of_triangles = no_triangles/0.7
964
965    log.critical('Estimated total number of triangles: %d'
966                 % total_number_of_triangles)
967    log.critical('Note: This is generally about 20%% '
968                 'less than the final amount')
969
970    return int(total_number_of_triangles)
971
972
973def decimate_polygon(polygon, factor=10):
974    """Reduce number of points in polygon by the specified
975    factor (default=10, hence the name of the function) such that
976    the extrema in both axes are preserved.
977
978##
979# @brief Reduce number of points in polygon by the specified factor.
980# @param polygon The polygon to reduce.
981# @param factor The factor to reduce polygon points by (default 10).
982# @note The extrema of both axes are preserved.
983
984    Return reduced polygon
985    """
986
987    # FIXME(Ole): This doesn't work at present,
988    # but it isn't critical either
989
990    # Find outer extent of polygon
991    num_polygon = ensure_numeric(polygon)
992    max_x = max(num_polygon[:, 0])
993    max_y = max(num_polygon[:, 1])
994    min_x = min(num_polygon[:, 0])
995    min_y = min(num_polygon[:, 1])
996
997    # Keep only some points making sure extrema are kept
998    reduced_polygon = []
999    for i, point in enumerate(polygon):
1000        if point in [min_x, max_x] and point in [min_y, max_y]:
1001            # Keep
1002            reduced_polygon.append(point)
1003        else:
1004            if len(reduced_polygon)*factor < i:
1005                reduced_polygon.append(point)
1006
1007    return reduced_polygon
1008
1009
1010def interpolate_polyline(data,
1011                         polyline_nodes,
1012                         gauge_neighbour_id,
1013                         interpolation_points=None,
1014                         rtol=1.0e-6,
1015                         atol=1.0e-8):
1016    """Interpolate linearly between values data on polyline nodes
1017    of a polyline to list of interpolation points.
1018
1019    data is the data on the polyline nodes.
1020
1021    Inputs:
1022      data: Vector or array of data at the polyline nodes.
1023      polyline_nodes: Location of nodes where data is available.
1024      gauge_neighbour_id: ?
1025      interpolation_points: Interpolate polyline data to these positions.
1026          List of coordinate pairs [x, y] of
1027          data points or an nx2 numeric array or a Geospatial_data object
1028      rtol, atol: Used to determine whether a point is on the polyline or not.
1029                  See point_on_line.
1030
1031    Output:
1032      Interpolated values at interpolation points
1033    """
1034
1035    if isinstance(interpolation_points, Geospatial_data):
1036        interpolation_points = interpolation_points.\
1037                                    get_data_points(absolute=True)
1038
1039    interpolated_values = num.zeros(len(interpolation_points), num.float)
1040
1041    data = ensure_numeric(data, num.float)
1042    polyline_nodes = ensure_numeric(polyline_nodes, num.float)
1043    interpolation_points = ensure_numeric(interpolation_points, num.float)
1044    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
1045
1046    num_nodes = polyline_nodes.shape    # Number of nodes in polyline
1047
1048    # Input sanity check
1049    assert_msg = 'interpolation_points are not given (interpolate.py)'
1050    assert interpolation_points is not None, assert_msg
1051
1052    assert_msg = 'function value must be specified at every interpolation node'
1053    assert data.shape == polyline_nodes.shape, assert_msg
1054
1055    assert_msg = 'Must define function value at one or more nodes'
1056    assert data.shape > 0, assert_msg
1057
1058    if num_nodes == 1:
1059        assert_msg = 'Polyline contained only one point. I need more. '
1060        assert_msg += str(data)
1061        raise Exception, assert_msg
1062    elif num_nodes > 1:
1063        _interpolate_polyline(data,
1064                              polyline_nodes,
1065                              gauge_neighbour_id,
1066                              interpolation_points,
1067                              interpolated_values,
1068                              rtol,
1069                              atol)
1070
1071
1072    return interpolated_values
1073
1074
1075def polylist2points_verts(polylist):
1076    """ Convert a list of polygons to discrete points and vertices.
1077    """
1078
1079    offset = 0
1080    points = []
1081    vertices = []
1082    for poly in polylist:
1083        points.extend(poly)
1084        vertices.extend([[i, i+1] for i in range(offset, offset+len(poly)-1)])
1085        offset += len(poly)
1086
1087    return points, vertices
1088
1089
1090################################################################################
1091# Initialise module
1092################################################################################
1093
1094from anuga.utilities import compile
1095if compile.can_use_C_extension('polygon_ext.c'):
1096    # Underlying C implementations can be accessed
1097    from polygon_ext import _point_on_line
1098    from polygon_ext import _separate_points_by_polygon
1099    from polygon_ext import _interpolate_polyline
1100    from polygon_ext import _is_inside_triangle
1101    #from polygon_ext import _intersection
1102
1103else:
1104    ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__
1105    ERROR_MSG += 'Make sure compile_all.py has been run as described in '
1106    ERROR_MSG += 'the ANUGA installation guide.'
1107    raise Exception(ERROR_MSG)
1108
1109
1110if __name__ == "__main__":
1111    pass
Note: See TracBrowser for help on using the repository browser.