source: anuga_core/source/anuga/utilities/polygon.py @ 6535

Last change on this file since 6535 was 6535, checked in by ole, 15 years ago

Wrote is_inside_triangle in C and tested it some more.

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