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

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

Implemented fix to ticket:314.

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