source: branches/numpy/anuga/utilities/polygon.py @ 6883

Last change on this file since 6883 was 6689, checked in by rwilson, 16 years ago

Back-merge from Numeric trunk to numpy branch.

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