# source:anuga_core/source/anuga/utilities/polygon.py@5570

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

Added input tests regarding polygons and points

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