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

Last change on this file since 5372 was 5372, checked in by jakeman, 16 years ago

updated mux2 boundary interpolation methods

File size: 29.0 KB
Line 
1#!/usr/bin/env python
2"""Polygon manipulations
3
4"""
5
6
7#try:
8#    from scipy import Float, Int, zeros, ones, array, concatenate, reshape, dot
9#except:
10#    #print 'Could not find scipy - using Numeric'
11
12from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose
13
14
15from math import sqrt
16from anuga.utilities.numerical_tools import ensure_numeric
17from anuga.geospatial_data.geospatial_data import ensure_absolute
18
19def point_on_line_py(point, line):
20    from Numeric import fabs
21    point = ensure_numeric(point)
22    line = ensure_numeric(line)
23
24
25    x=point[0];y=point[1]
26    x0=line[0,0];y0=line[0,1]
27    x1=line[1,0];y1=line[1,1]
28    #from pylab import plot,show
29    #plot(line[:,0],line[:,1])
30    #plot([x],[y],'o')
31    #show()
32    a0 = x - x0;
33    a1 = y - y0;
34
35    a_normal0 = a1;
36    a_normal1 = -a0;
37
38    b0 = x1 - x0;
39    b1 = y1 - y0;
40
41    #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ):
42    eps=200
43    #print 'normal',a_normal0*b0 + a_normal1*b1
44    if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ):
45        #for some reason (perhaps catastrophic cancellation) urs_boundary.py
46        #example only works for eps=2
47        #point is somewhere on the infinite extension of the line
48        # FIXME (Ole): Perhaps add a tolerance here instead of 0.0
49
50        len_a = sqrt(a0*a0 + a1*a1);
51        len_b = sqrt(b0*b0 + b1*b1);
52
53        if (a0*b0 + a1*b1 >= 0 and len_a <= len_b):
54            return 1
55        else:
56            return 0
57    else:
58        return 0
59
60
61def point_on_line(point, line):
62    """Determine whether a point is on a line segment
63
64    Input:
65        point is given by [x, y]
66        line is given by [x0, y0], [x1, y1]] or
67        the equivalent 2x2 Numeric array with each row corresponding to a point.
68
69    Output:
70       
71    """
72
73    point = ensure_numeric(point)
74    line = ensure_numeric(line)
75
76
77    res = _point_on_line(point[0], point[1],
78                         line[0,0], line[0,1],
79                         line[1,0], line[1,1])
80   
81    return bool(res)
82
83
84
85
86
87def intersection(line0, line1):
88    """Returns intersecting point between two line segments or None
89    (if parallel or no intersection is found).
90
91    However, if parallel lines coincide partly (i.e. shara a common segment,
92    the line segment where lines coincide is returned
93   
94
95    Inputs:
96        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
97                      A line can also be a 2x2 numeric array with each row
98                      corresponding to a point.
99
100
101    Output:
102        status, value
103
104        where status is interpreted as follows
105       
106        status == 0: no intersection with value set to None
107        status == 1: One intersection point found and returned in value as [x,y]
108        status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]]
109        status == 3: Lines would coincide but only if extended. Value set to None
110        status == 4: Lines are parallel with a fixed distance apart. Value set to None.
111   
112    """
113
114    # FIXME (Ole): Write this in C
115
116    line0 = ensure_numeric(line0, Float)
117    line1 = ensure_numeric(line1, Float)   
118
119    x0 = line0[0,0]; y0 = line0[0,1]
120    x1 = line0[1,0]; y1 = line0[1,1]
121
122    x2 = line1[0,0]; y2 = line1[0,1]
123    x3 = line1[1,0]; y3 = line1[1,1]
124
125    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
126    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
127    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
128       
129    if allclose(denom, 0.0):
130        # Lines are parallel - check if they coincide on a shared a segment
131
132        if allclose( [u0, u1], 0.0 ):
133            # We now know that the lines if continued coincide
134            # The remaining check will establish if the finite lines share a segment
135
136            line0_starts_on_line1 = line0_ends_on_line1 =\
137            line1_starts_on_line0 = line1_ends_on_line0 = False
138               
139            if point_on_line([x0, y0], line1):
140                line0_starts_on_line1 = True
141
142            if point_on_line([x1, y1], line1):
143                line0_ends_on_line1 = True
144 
145            if point_on_line([x2, y2], line0):
146                line1_starts_on_line0 = True
147
148            if point_on_line([x3, y3], line0):
149                line1_ends_on_line0 = True                               
150
151            if not(line0_starts_on_line1 or line0_ends_on_line1\
152               or line1_starts_on_line0 or line1_ends_on_line0):
153                # Lines are parallel and would coincide if extended, but not as they are.
154                return 3, None
155
156
157            # One line fully included in the other. Use direction of included line
158            if line0_starts_on_line1 and line0_ends_on_line1:
159                # Shared segment is line0 fully included in line1
160                segment = array([[x0, y0], [x1, y1]])               
161
162            if line1_starts_on_line0 and line1_ends_on_line0:
163                # Shared segment is line1 fully included in line0
164                segment = array([[x2, y2], [x3, y3]])
165           
166
167            # Overlap with lines are oriented the same way
168            if line0_starts_on_line1 and line1_ends_on_line0:
169                # Shared segment from line0 start to line 1 end
170                segment = array([[x0, y0], [x3, y3]])
171
172            if line1_starts_on_line0 and line0_ends_on_line1:
173                # Shared segment from line1 start to line 0 end
174                segment = array([[x2, y2], [x1, y1]])                               
175
176
177            # Overlap in opposite directions - use direction of line0
178            if line0_starts_on_line1 and line1_starts_on_line0:
179                # Shared segment from line0 start to line 1 end
180                segment = array([[x0, y0], [x2, y2]])
181
182            if line0_ends_on_line1 and line1_ends_on_line0:
183                # Shared segment from line0 start to line 1 end
184                segment = array([[x3, y3], [x1, y1]])               
185
186               
187            return 2, segment
188        else:
189            # Lines are parallel but they don't coincide
190            return 4, None #FIXME (Ole): Add distance here instead of None
191           
192    else:
193        # Lines are not parallel or coinciding
194        u0 = u0/denom
195        u1 = u1/denom       
196
197        x = x0 + u0*(x1-x0)
198        y = y0 + u0*(y1-y0)
199
200        # Sanity check - can be removed to speed up if needed
201        assert allclose(x, x2 + u1*(x3-x2))
202        assert allclose(y, y2 + u1*(y3-y2))       
203
204        # Check if point found lies within given line segments
205        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 
206            # We have intersection
207
208            return 1, array([x, y])
209        else:
210            # No intersection
211            return 0, None
212
213
214def NEW_C_intersection(line0, line1):
215    #FIXME(Ole): To write in C
216    """Returns intersecting point between two line segments or None
217    (if parallel or no intersection is found).
218
219    However, if parallel lines coincide partly (i.e. shara a common segment,
220    the line segment where lines coincide is returned
221   
222
223    Inputs:
224        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
225                      A line can also be a 2x2 numeric array with each row
226                      corresponding to a point.
227
228
229    Output:
230        status, value
231
232        where status is interpreted as follows
233       
234        status == 0: no intersection with value set to None
235        status == 1: One intersection point found and returned in value as [x,y]
236        status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]]
237        status == 3: Lines would coincide but only if extended. Value set to None
238        status == 4: Lines are parallel with a fixed distance apart. Value set to None.
239   
240    """
241
242
243    line0 = ensure_numeric(line0, Float)
244    line1 = ensure_numeric(line1, Float)   
245
246    status, value = _intersection(line0[0,0], line0[0,1],
247                                  line0[1,0], line0[1,1],
248                                  line1[0,0], line1[0,1],
249                                  line1[1,0], line1[1,1])
250
251    return status, value
252
253
254
255
256def is_inside_polygon(point, polygon, closed=True, verbose=False):
257    """Determine if one point is inside a polygon
258
259    See inside_polygon for more details
260    """
261
262    indices = inside_polygon(point, polygon, closed, verbose)
263
264    if indices.shape[0] == 1:
265        return True
266    elif indices.shape[0] == 0:
267        return False
268    else:
269        msg = 'is_inside_polygon must be invoked with one point only'
270        raise msg
271   
272
273def inside_polygon(points, polygon, closed=True, verbose=False):
274    """Determine points inside a polygon
275
276       Functions inside_polygon and outside_polygon have been defined in
277       terms af separate_by_polygon which will put all inside indices in
278       the first part of the indices array and outside indices in the last
279
280       See separate_points_by_polygon for documentation
281
282       points and polygon can be a geospatial instance,
283       a list or a numeric array
284    """
285
286    #if verbose: print 'Checking input to inside_polygon'
287
288    try:
289        points = ensure_absolute(points)
290    except NameError, e:
291        raise NameError, e
292    except:
293        # If this fails it is going to be because the points can't be
294        # converted to a numeric array.
295        msg = 'Points could not be converted to Numeric array' 
296        raise msg
297
298    try:
299        polygon = ensure_absolute(polygon)
300    except NameError, e:
301        raise NameError, e
302    except:
303        # If this fails it is going to be because the points can't be
304        # converted to a numeric array.
305        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
306        raise msg
307
308    if len(points.shape) == 1:
309        # Only one point was passed in. Convert to array of points
310        points = reshape(points, (1,2))
311
312    indices, count = separate_points_by_polygon(points, polygon,
313                                                closed=closed,
314                                                verbose=verbose)
315
316    # Return indices of points inside polygon
317    return indices[:count]
318
319
320
321def is_outside_polygon(point, polygon, closed=True, verbose=False,
322                       points_geo_ref=None, polygon_geo_ref=None):
323    """Determine if one point is outside a polygon
324
325    See outside_polygon for more details
326    """
327
328    indices = outside_polygon(point, polygon, closed, verbose)
329                              #points_geo_ref, polygon_geo_ref)
330
331    if indices.shape[0] == 1:
332        return True
333    elif indices.shape[0] == 0:
334        return False
335    else:
336        msg = 'is_outside_polygon must be invoked with one point only'
337        raise msg
338   
339
340def outside_polygon(points, polygon, closed = True, verbose = False):
341    """Determine points outside a polygon
342
343       Functions inside_polygon and outside_polygon have been defined in
344       terms af separate_by_polygon which will put all inside indices in
345       the first part of the indices array and outside indices in the last
346
347       See separate_points_by_polygon for documentation
348    """
349
350    #if verbose: print 'Checking input to outside_polygon'
351    try:
352        points = ensure_numeric(points, Float)
353    except NameError, e:
354        raise NameError, e
355    except:
356        msg = 'Points could not be converted to Numeric array'
357        raise msg
358
359    try:
360        polygon = ensure_numeric(polygon, Float)
361    except NameError, e:
362        raise NameError, e
363    except:
364        msg = 'Polygon could not be converted to Numeric array'
365        raise msg
366
367
368    if len(points.shape) == 1:
369        # Only one point was passed in. Convert to array of points
370        points = reshape(points, (1,2))
371
372    indices, count = separate_points_by_polygon(points, polygon,
373                                                closed=closed,
374                                                verbose=verbose)
375
376    # Return indices of points outside polygon
377    if count == len(indices):
378        # No points are outside
379        return array([])
380    else:
381        return indices[count:][::-1]  #return reversed
382       
383
384def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
385    """Determine points inside and outside a polygon
386
387       See separate_points_by_polygon for documentation
388
389       Returns an array of points inside and an array of points outside the polygon
390    """
391
392    #if verbose: print 'Checking input to outside_polygon'
393    try:
394        points = ensure_numeric(points, Float)
395    except NameError, e:
396        raise NameError, e
397    except:
398        msg = 'Points could not be converted to Numeric array'
399        raise msg
400
401    try:
402        polygon = ensure_numeric(polygon, Float)
403    except NameError, e:
404        raise NameError, e
405    except:
406        msg = 'Polygon could not be converted to Numeric array'
407        raise msg
408
409    if len(points.shape) == 1:
410        # Only one point was passed in. Convert to array of points
411        points = reshape(points, (1,2))
412
413
414    indices, count = separate_points_by_polygon(points, polygon,
415                                                closed=closed,
416                                                verbose=verbose)
417   
418    # Returns indices of points inside and indices of points outside
419    # the polygon
420
421    if count == len(indices):
422        # No points are outside
423        return indices[:count],[]
424    else:
425        return  indices[:count], indices[count:][::-1]  #return reversed
426
427
428def separate_points_by_polygon(points, polygon,
429                               closed = True, verbose = False):
430    """Determine whether points are inside or outside a polygon
431
432    Input:
433       points - Tuple of (x, y) coordinates, or list of tuples
434       polygon - list of vertices of polygon
435       closed - (optional) determine whether points on boundary should be
436       regarded as belonging to the polygon (closed = True)
437       or not (closed = False)
438
439    Outputs:
440       indices: array of same length as points with indices of points falling
441       inside the polygon listed from the beginning and indices of points
442       falling outside listed from the end.
443
444       count: count of points falling inside the polygon
445
446       The indices of points inside are obtained as indices[:count]
447       The indices of points outside are obtained as indices[count:]
448
449
450    Examples:
451       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
452
453       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
454       will return the indices [0, 2, 1] and count == 2 as only the first
455       and the last point are inside the unit square
456
457    Remarks:
458       The vertices may be listed clockwise or counterclockwise and
459       the first point may optionally be repeated.
460       Polygons do not need to be convex.
461       Polygons can have holes in them and points inside a hole is
462       regarded as being outside the polygon.
463
464    Algorithm is based on work by Darel Finley,
465    http://www.alienryderflex.com/polygon/
466
467    Uses underlying C-implementation in polygon_ext.c
468    """
469
470
471    #if verbose: print 'Checking input to separate_points_by_polygon'
472
473
474    #Input checks
475
476    assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
477    assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
478
479
480    try:
481        points = ensure_numeric(points, Float)
482    except NameError, e:
483        raise NameError, e
484    except:
485        msg = 'Points could not be converted to Numeric array'
486        raise msg
487
488    #if verbose: print 'Checking input to separate_points_by_polygon 2'
489    try:
490        polygon = ensure_numeric(polygon, Float)
491    except NameError, e:
492        raise NameError, e
493    except:
494        msg = 'Polygon could not be converted to Numeric array'
495        raise msg
496
497    #if verbose: print 'check'
498
499    assert len(polygon.shape) == 2,\
500       'Polygon array must be a 2d array of vertices'
501
502    assert polygon.shape[1] == 2,\
503       'Polygon array must have two columns'
504
505    assert len(points.shape) == 2,\
506       'Points array must be a 2d array'
507
508    assert points.shape[1] == 2,\
509       'Points array must have two columns'
510
511    N = polygon.shape[0] #Number of vertices in polygon
512    M = points.shape[0]  #Number of points
513
514
515    indices = zeros( M, Int )
516
517    count = _separate_points_by_polygon(points, polygon, indices,
518                                        int(closed), int(verbose))
519
520    if verbose: print 'Found %d points (out of %d) inside polygon'\
521       %(count, M)
522    return indices, count
523
524
525def polygon_area(polygon):
526    """ Determin area of arbitrary polygon
527    Reference
528    http://mathworld.wolfram.com/PolygonArea.html
529    """
530   
531    n = len(polygon)
532    poly_area = 0.0
533
534    for i in range(n):
535        pti = polygon[i]
536        if i == n-1:
537            pt1 = polygon[0]
538        else:
539            pt1 = polygon[i+1]
540        xi = pti[0]
541        yi1 = pt1[1]
542        xi1 = pt1[0]
543        yi = pti[1]
544        poly_area += xi*yi1 - xi1*yi
545       
546    return abs(poly_area/2)
547
548def plot_polygons(polygons_points, style=None, 
549                  figname=None, label=None, verbose=False):
550   
551    """ Take list of polygons and plot.
552
553    Inputs:
554
555    polygons         - list of polygons
556
557    style            - style list corresponding to each polygon
558                     - for a polygon, use 'line'
559                     - for points falling outside a polygon, use 'outside'
560                       
561    figname          - name to save figure to
562
563    label            - title for plot
564
565    Outputs:
566
567    - list of min and max of x and y coordinates
568    - plot of polygons
569    """ 
570
571    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close, title
572
573    assert type(polygons_points) == list,\
574               'input must be a list of polygons and/or points'
575               
576    ion()
577    hold(True)
578
579    minx = 1e10
580    maxx = 0.0
581    miny = 1e10
582    maxy = 0.0
583
584    if label is None: label = ''
585
586    n = len(polygons_points)
587    colour = []
588    if style is None:
589        style_type = 'line' 
590        style = []
591        for i in range(n):
592            style.append(style_type)
593            colour.append('b-')
594    else:
595        for s in style:
596            if s == 'line': colour.append('b-')           
597            if s == 'outside': colour.append('r.')
598            if s <> 'line':
599                if s <> 'outside':
600                    colour.append('g.')
601           
602    for i, item in enumerate(polygons_points):
603        x, y = poly_xy(item)
604        if min(x) < minx: minx = min(x)
605        if max(x) > maxx: maxx = max(x)
606        if min(y) < miny: miny = min(y)
607        if max(y) > maxy: maxy = max(y)
608        plot(x,y,colour[i])
609        xlabel('x')
610        ylabel('y')
611        title(label)
612
613    #raw_input('wait 1')
614    #FIXME(Ole): This makes for some strange scalings sometimes.
615    #if minx <> 0:
616    #    axis([minx*0.9,maxx*1.1,miny*0.9,maxy*1.1])
617    #else:
618    #    if miny == 0:
619    #        axis([-maxx*.01,maxx*1.1,-maxy*0.01,maxy*1.1])
620    #    else:
621    #        axis([-maxx*.01,maxx*1.1,miny*0.9,maxy*1.1])
622
623    if figname is not None:
624        savefig(figname)
625    else:
626        savefig('test_image')
627
628    close('all')
629
630    vec = [minx,maxx,miny,maxy]
631
632    return vec
633
634def poly_xy(polygon, verbose=False):
635    """ this is used within plot_polygons so need to duplicate
636        the first point so can have closed polygon in plot
637    """
638
639    #if verbose: print 'Checking input to poly_xy'
640
641    try:
642        polygon = ensure_numeric(polygon, Float)
643    except NameError, e:
644        raise NameError, e
645    except:
646        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
647        raise msg
648
649    x = polygon[:,0]
650    y = polygon[:,1]
651    x = concatenate((x, [polygon[0,0]]), axis = 0)
652    y = concatenate((y, [polygon[0,1]]), axis = 0)
653   
654    return x, y
655   
656#    x = []
657#    y = []
658#    n = len(poly)
659#    firstpt = poly[0]
660#    for i in range(n):
661#        thispt = poly[i]
662#        x.append(thispt[0])
663#        y.append(thispt[1])
664
665#    x.append(firstpt[0])
666#    y.append(firstpt[1])
667   
668#    return x, y
669
670class Polygon_function:
671    """Create callable object f: x,y -> z, where a,y,z are vectors and
672    where f will return different values depending on whether x,y belongs
673    to specified polygons.
674
675    To instantiate:
676
677       Polygon_function(polygons)
678
679    where polygons is a list of tuples of the form
680
681      [ (P0, v0), (P1, v1), ...]
682
683      with Pi being lists of vertices defining polygons and vi either
684      constants or functions of x,y to be applied to points with the polygon.
685
686    The function takes an optional argument, default which is the value
687    (or function) to used for points not belonging to any polygon.
688    For example:
689
690       Polygon_function(polygons, default = 0.03)
691
692    If omitted the default value will be 0.0
693
694    Note: If two polygons overlap, the one last in the list takes precedence
695
696    Coordinates specified in the call are assumed to be relative to the
697    origin (georeference) e.g. used by domain.
698    By specifying the optional argument georeference,
699    all points are made relative.
700
701    FIXME: This should really work with geo_spatial point sets.
702    """
703
704    def __init__(self, regions, default=0.0, geo_reference=None):
705
706        try:
707            len(regions)
708        except:
709            msg = 'Polygon_function takes a list of pairs (polygon, value).'
710            msg += 'Got %s' %polygons
711            raise msg
712
713
714        T = regions[0]
715        try:
716            a = len(T)
717        except:
718            msg = 'Polygon_function takes a list of pairs (polygon, value).'
719            msg += 'Got %s' %polygons
720            raise msg
721
722        assert a == 2, 'Must have two component each: %s' %T
723
724
725        if geo_reference is None:
726            from anuga.coordinate_transforms.geo_reference import Geo_reference
727            geo_reference = Geo_reference()
728
729
730        self.default = default
731
732        # Make points in polygons relative to geo_reference
733        self.regions = []
734        for polygon, value in regions:
735            P = geo_reference.change_points_geo_ref(polygon)
736            self.regions.append( (P, value) )
737
738
739
740
741    def __call__(self, x, y):
742        x = array(x).astype(Float)
743        y = array(y).astype(Float)
744
745        N = len(x)
746        assert len(y) == N
747
748        points = concatenate( (reshape(x, (N, 1)),
749                               reshape(y, (N, 1))), axis=1 )
750
751        if callable(self.default):
752            z = self.default(x,y)
753        else:
754            z = ones(N, Float) * self.default
755
756        for polygon, value in self.regions:
757            indices = inside_polygon(points, polygon)
758
759            # FIXME: This needs to be vectorised
760            if callable(value):
761                for i in indices:
762                    xx = array([x[i]])
763                    yy = array([y[i]])
764                    z[i] = value(xx, yy)[0]
765            else:
766                for i in indices:
767                    z[i] = value
768
769        return z
770
771
772def read_polygon(filename, split=','):
773    """Read points assumed to form a polygon.
774       There must be exactly two numbers in each line separated by a comma.
775       No header.
776    """
777
778    #Get polygon
779    fid = open(filename)
780    lines = fid.readlines()
781    fid.close()
782    polygon = []
783    for line in lines:
784        fields = line.split(split)
785        polygon.append( [float(fields[0]), float(fields[1])] )
786
787    return polygon
788
789
790def write_polygon(polygon, filename=None):
791    """Write polygon to csv file.
792       There will be exactly two numbers, easting and northing,
793       in each line separated by a comma.
794       
795       No header.   
796    """
797
798    fid = open(filename, 'w')
799    for point in polygon:
800        fid.write('%f, %f\n' %point)
801    fid.close()
802   
803
804def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
805    """Populate given polygon with uniformly distributed points.
806
807    Input:
808       polygon - list of vertices of polygon
809       number_of_points - (optional) number of points
810       seed - seed for random number generator (default=None)
811       exclude - list of polygons (inside main polygon) from where points should be excluded
812
813    Output:
814       points - list of points inside polygon
815
816    Examples:
817       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
818       will return five randomly selected points inside the unit square
819    """
820
821    from random import uniform, seed as seed_function
822
823    seed_function(seed)
824
825    points = []
826
827    # Find outer extent of polygon
828    max_x = min_x = polygon[0][0]
829    max_y = min_y = polygon[0][1]
830    for point in polygon[1:]:
831        x = point[0]
832        if x > max_x: max_x = x
833        if x < min_x: min_x = x
834        y = point[1]
835        if y > max_y: max_y = y
836        if y < min_y: min_y = y
837
838
839    while len(points) < number_of_points:
840        x = uniform(min_x, max_x)
841        y = uniform(min_y, max_y)
842
843        append = False
844        if is_inside_polygon([x,y], polygon):
845
846            append = True
847
848            #Check exclusions
849            if exclude is not None:
850                for ex_poly in exclude:
851                    if is_inside_polygon([x,y], ex_poly):
852                        append = False
853
854
855        if append is True:
856            points.append([x,y])
857
858    return points
859
860
861def point_in_polygon(polygon, delta=1e-8):
862    """Return a point inside a given polygon which will be close to the
863    polygon edge.
864
865    Input:
866       polygon - list of vertices of polygon
867       delta - the square root of 2 * delta is the maximum distance from the
868       polygon points and the returned point.
869    Output:
870       points - a point inside polygon
871
872       searches in all diagonals and up and down (not left and right)
873    """
874    import exceptions
875    class Found(exceptions.Exception): pass
876
877    point_in = False
878    while not point_in:
879        try:
880            for poly_point in polygon: #[1:]:
881                for x_mult in range (-1,2):
882                    for y_mult in range (-1,2):
883                        x = poly_point[0]
884                        y = poly_point[1]
885                        if x == 0:
886                            x_delta = x_mult*delta
887                        else:
888                            x_delta = x+x_mult*x*delta
889
890                        if y == 0:
891                            y_delta = y_mult*delta
892                        else:
893                            y_delta = y+y_mult*y*delta
894
895                        point = [x_delta, y_delta]
896                        #print "point",point
897                        if is_inside_polygon(point, polygon, closed=False):
898                            raise Found
899        except Found:
900            point_in = True
901        else:
902            delta = delta*0.1
903    return point
904
905
906def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
907    """Calculate the approximate number of triangles inside the
908    bounding polygon and the other interior regions
909
910    Polygon areas are converted to square Kms
911
912    FIXME: Add tests for this function
913    """
914   
915    from anuga.utilities.polygon import polygon_area
916
917
918    # TO DO check if any of the regions fall inside one another
919
920    print '----------------------------------------------------------------------------'
921    print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
922    print '----------------------------------------------------------------------------'   
923       
924    no_triangles = 0.0
925    area = polygon_area(bounding_poly)
926   
927    for poly, resolution in interior_regions:
928        this_area = polygon_area(poly)
929        this_triangles = this_area/resolution
930        no_triangles += this_triangles
931        area -= this_area
932       
933        print 'Interior ',
934        print ('%.0f' %resolution).ljust(25),
935        print ('%.2f' %(this_area/1000000)).ljust(19),
936        print '%d' %(this_triangles)
937       
938    bound_triangles = area/remainder_res
939    no_triangles += bound_triangles
940
941    print 'Bounding ',
942    print ('%.0f' %remainder_res).ljust(25),
943    print ('%.2f' %(area/1000000)).ljust(19),
944    print '%d' %(bound_triangles)   
945
946    total_number_of_triangles = no_triangles/0.7
947
948    print 'Estimated total number of triangles: %d' %total_number_of_triangles
949    print 'Note: This is generally about 20% less than the final amount'   
950
951    return int(total_number_of_triangles)
952
953
954def decimate_polygon(polygon, factor=10):
955    """Reduce number of points in polygon by the specified
956    factor (default=10, hence the name of the function) such that
957    the extrema in both axes are preserved.
958
959    Return reduced polygon
960    """
961
962    # FIXME(Ole): This doesn't work at present,
963    # but it isn't critical either
964
965    # Find outer extent of polygon
966    num_polygon = ensure_numeric(polygon)
967    max_x = max(num_polygon[:,0])
968    max_y = max(num_polygon[:,1])
969    min_x = min(num_polygon[:,0])
970    min_y = min(num_polygon[:,1])       
971
972    # Keep only some points making sure extrema are kept
973    reduced_polygon = []   
974    for i, point in enumerate(polygon):
975        x = point[0]
976        y = point[1]       
977        if x in [min_x, max_x] and y in [min_y, max_y]:
978            # Keep
979            reduced_polygon.append(point)
980        else:
981            if len(reduced_polygon)*factor < i:
982                reduced_polygon.append(point)               
983
984    return reduced_polygon
985
986##############################################
987#Initialise module
988
989from anuga.utilities.compile import can_use_C_extension
990if can_use_C_extension('polygon_ext.c'):
991    # Underlying C implementations can be accessed
992    from polygon_ext import _point_on_line
993    from polygon_ext import _separate_points_by_polygon
994    #from polygon_ext import _intersection
995
996else:
997    msg = 'C implementations could not be accessed by %s.\n ' %__file__
998    msg += 'Make sure compile_all.py has been run as described in '
999    msg += 'the ANUGA installation guide.'
1000    raise Exception, msg
1001
1002
1003if __name__ == "__main__":
1004    pass
Note: See TracBrowser for help on using the repository browser.