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

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

Added new function is_inside_triangle. This helps the fitting functions.

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