Changeset 1911

Timestamp:

Oct 14, 2005, 6:44:45 AM (17 years ago)

Author:

ole

Message:

Refactored pyvolution to use polygon functionality from new utilities package

Location:

inundation

Files:

• pyvolution/combine_pts.py (modified) (1 diff)
• pyvolution/data_manager.py (modified) (6 diffs)
• pyvolution/least_squares.py (modified) (2 diffs)
• pyvolution/test_mesh.py (modified) (5 diffs)
• pyvolution/test_util.py (modified) (2 diffs)
• pyvolution/util.py (modified) (4 diffs)
• utilities/polygon.py (modified) (2 diffs)

inundation/pyvolution/combine_pts.py

@@ -8 +8 @@
    Geoscience Australia, 2005.   
 """
-from util import outside_polygon, inside_polygon
+from utilities.polygon import outside_polygon, inside_polygon
 from Numeric import take, concatenate
 import time

inundation/pyvolution/data_manager.py

@@ -1261 +1261 @@
 
     from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 
-
+    from utilities.polygon import inside_polygon
+    
     msg = 'Format must be either asc or ers'
     assert format.lower() in ['asc', 'ers'], msg
@@ -1399 +1400 @@
     #Interpolate
     from least_squares import Interpolation
-    from util import inside_polygon
+    
 
     #FIXME: This should be done with precrop = True (?), otherwise it'll
@@ -2880 +2881 @@
     from Numeric import array, Float, concatenate, NewAxis, zeros,\
          sometrue
-
+    from utilities.polygon import inside_polygon
 
     #FIXME: Should be variable
@@ -2984 +2985 @@
     #Interpolate
     from least_squares import Interpolation
-    from util import inside_polygon
+    
 
     #FIXME: This should be done with precrop = True, otherwise it'll
@@ -3097 +3098 @@
     from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 
     import ermapper_grids
+    from utilities.polygon import inside_polygon    
 
     header = {}
@@ -3202 +3204 @@
     #Interpolate
     from least_squares import Interpolation
-    from util import inside_polygon
+    
 
     #FIXME: This should be done with precrop = True (?), otherwise it'll

inundation/pyvolution/least_squares.py

@@ -434 +434 @@
 
         from quad import build_quadtree
-        from util import ensure_numeric
+        from utilities.numerical_tools import ensure_numeric
+        from utilities.polygon import inside_polygon
+        
 
         if data_origin is None:
@@ -489 +491 @@
         if precrop is True:
             from Numeric import take
-            from util import inside_polygon
 
             if verbose: print 'Getting boundary polygon'

inundation/pyvolution/test_mesh.py

@@ -10 +10 @@
 from config import epsilon
 from Numeric import allclose, array
+
+from utilities.polygon import inside_polygon
 
 def distance(x, y):
@@ -705 +707 @@
         from mesh import Mesh
         from Numeric import zeros, Float
-        from util import inside_polygon
 
         #Create basic mesh
@@ -726 +727 @@
         from mesh import Mesh
         from Numeric import zeros, Float
-        from util import inside_polygon
+        
 
         #Points
@@ -766 +767 @@
         from mesh import Mesh
         from Numeric import zeros, Float
-        from util import inside_polygon
+
 
         #Points
@@ -807 +808 @@
         from mesh import Mesh
         from Numeric import zeros, Float
-        from util import inside_polygon
         from mesh_factory import rectangular        
 

inundation/pyvolution/test_util.py

@@ -172 +172 @@
         assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
 
-    def test_ensure_numeric(self):
-        from util import ensure_numeric
-        from Numeric import ArrayType, Float, array
-
-        A = [1,2,3,4]
-        B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'l'
-        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
-
-
-        A = [1,2,3.14,4]
-        B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
-        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
-
-
-        A = [1,2,3,4]
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
-        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
-
-
-        A = [1,2,3,4]
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'
-        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
-
-
-        A = array([1,2,3,4])
-        B = ensure_numeric(A)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'l'        
-        assert A == B    
-        assert A is B   #Same object
-
-
-        A = array([1,2,3,4])
-        B = ensure_numeric(A, Float)
-        assert type(B) == ArrayType
-        assert B.typecode() == 'd'        
-        assert A == B    
-        assert A is not B   #Not the same object        
-
 
 
@@ -1020 +973 @@
 
 
-    #Polygon stuff
-    def test_polygon_function_constants(self):
-        p1 = [[0,0], [10,0], [10,10], [0,10]]
-        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-
-        f = Polygon_function( [(p1, 1.0)] )
-        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
-        assert allclose(z, [1,1,0,0])
-
-
-        f = Polygon_function( [(p2, 2.0)] )
-        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
-        assert allclose(z, [2,0,0,2])
-
-
-        #Combined
-        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
-        z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,1,0,2])
-
-
-    def test_polygon_function_callable(self):
-        """Check that values passed into Polygon_function can be callable
-        themselves.
-        """
-        p1 = [[0,0], [10,0], [10,10], [0,10]]
-        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-
-        f = Polygon_function( [(p1, test_function)] )
-        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
-        assert allclose(z, [10,14,0,0])
-
-        #Combined
-        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
-        z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,14,0,2])
-
-
-        #Combined w default
-        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
-        z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,14,3.14,2])
-
-
-        #Combined w default func
-        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
-                              default = test_function)
-        z = f([5, 5, 27, 35], [5, 9, 8, -5])
-        assert allclose(z, [2,14,35,2])
-
-
-    def test_point_on_line(self):
-
-        #Endpoints first
-        assert point_on_line( 0, 0, 0,0, 1,0 )
-        assert point_on_line( 1, 0, 0,0, 1,0 )
-
-        #Then points on line
-        assert point_on_line( 0.5, 0, 0,0, 1,0 )
-        assert point_on_line( 0, 0.5, 0,1, 0,0 )
-        assert point_on_line( 1, 0.5, 1,1, 1,0 )
-        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
-
-        #Then points not on line
-        assert not point_on_line( 0.5, 0, 0,1, 1,1 )
-        assert not point_on_line( 0, 0.5, 0,0, 1,1 )
-
-        #From real example that failed
-        assert not point_on_line( 40,50, 40,20, 40,40 )
-
-
-        #From real example that failed
-        assert not point_on_line( 40,19, 40,20, 40,40 )
-
-
-
-
-    def test_inside_polygon_main(self):
-
-
-        #Simplest case: Polygon is the unit square
-        polygon = [[0,0], [1,0], [1,1], [0,1]]
-
-        assert inside_polygon( (0.5, 0.5), polygon )
-        assert not inside_polygon( (0.5, 1.5), polygon )
-        assert not inside_polygon( (0.5, -0.5), polygon )
-        assert not inside_polygon( (-0.5, 0.5), polygon )
-        assert not inside_polygon( (1.5, 0.5), polygon )
-
-        #Try point on borders
-        assert inside_polygon( (1., 0.5), polygon, closed=True)
-        assert inside_polygon( (0.5, 1), polygon, closed=True)
-        assert inside_polygon( (0., 0.5), polygon, closed=True)
-        assert inside_polygon( (0.5, 0.), polygon, closed=True)
-
-        assert not inside_polygon( (0.5, 1), polygon, closed=False)
-        assert not inside_polygon( (0., 0.5), polygon, closed=False)
-        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
-        assert not inside_polygon( (1., 0.5), polygon, closed=False)
-
-
-
-        #From real example (that failed)
-        polygon = [[20,20], [40,20], [40,40], [20,40]]
-        points = [ [40, 50] ]
-        res = inside_polygon(points, polygon)
-        assert len(res) == 0
-
-        polygon = [[20,20], [40,20], [40,40], [20,40]]
-        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
-        res = inside_polygon(points, polygon)
-        assert len(res) == 2
-        assert allclose(res, [0,1])
-
-
-
-        #More convoluted and non convex polygon
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        assert inside_polygon( (0.5, 0.5), polygon )
-        assert inside_polygon( (1, -0.5), polygon )
-        assert inside_polygon( (1.5, 0), polygon )
-
-        assert not inside_polygon( (0.5, 1.5), polygon )
-        assert not inside_polygon( (0.5, -0.5), polygon )
-
-
-        #Very convoluted polygon
-        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-        assert inside_polygon( (5, 5), polygon )
-        assert inside_polygon( (17, 7), polygon )
-        assert inside_polygon( (27, 2), polygon )
-        assert inside_polygon( (35, -5), polygon )
-        assert not inside_polygon( (15, 7), polygon )
-        assert not inside_polygon( (24, 3), polygon )
-        assert not inside_polygon( (25, -10), polygon )
-
-
-
-        #Another combination (that failed)
-        polygon = [[0,0], [10,0], [10,10], [0,10]]
-        assert inside_polygon( (5, 5), polygon )
-        assert inside_polygon( (7, 7), polygon )
-        assert not inside_polygon( (-17, 7), polygon )
-        assert not inside_polygon( (7, 17), polygon )
-        assert not inside_polygon( (17, 7), polygon )
-        assert not inside_polygon( (27, 8), polygon )
-        assert not inside_polygon( (35, -5), polygon )
-
-
-
-
-        #Multiple polygons
-
-        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
-                   [10,10], [11,10], [11,11], [10,11], [10,10]]
-        assert inside_polygon( (0.5, 0.5), polygon )
-        assert inside_polygon( (10.5, 10.5), polygon )
-
-        #FIXME: Fails if point is 5.5, 5.5
-        assert not inside_polygon( (0, 5.5), polygon )
-
-        #Polygon with a hole
-        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
-                   [0,0], [1,0], [1,1], [0,1], [0,0]]
-
-        assert inside_polygon( (0, -0.5), polygon )
-        assert not inside_polygon( (0.5, 0.5), polygon )
-
-    def test_inside_polygon_vector_version(self):
-        #Now try the vector formulation returning indices
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        res = inside_polygon( points, polygon, verbose=False )
-
-        assert allclose( res, [0,1,2] )
-
-    def test_outside_polygon(self):
-        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
-
-        assert not outside_polygon( [0.5, 0.5], U )
-        #evaluate to False as the point 0.5, 0.5 is inside the unit square
-        
-        assert outside_polygon( [1.5, 0.5], U )
-        #evaluate to True as the point 1.5, 0.5 is outside the unit square
-        
-        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-        assert allclose( indices, [1] )
-        
-        #One more test of vector formulation returning indices
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        res = outside_polygon( points, polygon )
-
-        assert allclose( res, [3, 4] )
-
-
-
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        res = outside_polygon( points, polygon )
-
-        assert allclose( res, [0, 4, 5] )        
-     
-    def test_outside_polygon2(self):
-        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
-   
-        assert not outside_polygon( [0.5, 1.0], U, closed = True )
-        #evaluate to False as the point 0.5, 1.0 is inside the unit square
-        
-        assert outside_polygon( [0.5, 1.0], U, closed = False )
-        #evaluate to True as the point 0.5, 1.0 is outside the unit square
-
-    def test_separate_points_by_polygon(self):
-        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
-
-        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-        assert allclose( indices, [0,2,1] )
-        assert count == 2
-        
-        #One more test of vector formulation returning indices
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        res, count = separate_points_by_polygon( points, polygon )
-
-        assert allclose( res, [0,1,2,4,3] )
-        assert count == 3
-
-
-        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        res, count = separate_points_by_polygon( points, polygon )
-
-        assert allclose( res, [1,2,3,5,4,0] )        
-        assert count == 3
-        
-
-    def test_populate_polygon(self):
-
-        polygon = [[0,0], [1,0], [1,1], [0,1]]
-        points = populate_polygon(polygon, 5)
-
-        assert len(points) == 5
-        for point in points:
-            assert inside_polygon(point, polygon)
-
-
-        #Very convoluted polygon
-        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-
-        points = populate_polygon(polygon, 5)
-
-        assert len(points) == 5
-        for point in points:
-            assert inside_polygon(point, polygon)
-
-
-    def test_populate_polygon_with_exclude(self):
-        
-
-        polygon = [[0,0], [1,0], [1,1], [0,1]]
-        ex_poly = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
-        points = populate_polygon(polygon, 5, exclude = [ex_poly])
-
-        assert len(points) == 5
-        for point in points:
-            assert inside_polygon(point, polygon)
-            assert not inside_polygon(point, ex_poly)            
-
-
-        #overlap
-        polygon = [[0,0], [1,0], [1,1], [0,1]]
-        ex_poly = [[-1,-1], [0.5,0], [0.5, 0.5], [-1,0.5]]
-        points = populate_polygon(polygon, 5, exclude = [ex_poly])
-
-        assert len(points) == 5
-        for point in points:
-            assert inside_polygon(point, polygon)
-            assert not inside_polygon(point, ex_poly)                        
-        
-        #Multiple
-        polygon = [[0,0], [1,0], [1,1], [0,1]]
-        ex_poly1 = [[0,0], [0.5,0], [0.5, 0.5], [0,0.5]] #SW quarter
-        ex_poly2 = [[0.5,0.5], [0.5,1], [1, 1], [1,0.5]] #NE quarter        
-        
-        points = populate_polygon(polygon, 20, exclude = [ex_poly1, ex_poly2])
-
-        assert len(points) == 20
-        for point in points:
-            assert inside_polygon(point, polygon)
-            assert not inside_polygon(point, ex_poly1)
-            assert not inside_polygon(point, ex_poly2)                                
-        
-
-        #Very convoluted polygon
-        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
-        ex_poly = [[-1,-1], [5,0], [5, 5], [-1,5]]
-        points = populate_polygon(polygon, 20, exclude = [ex_poly])
-        
-        assert len(points) == 20
-        for point in points:
-            assert inside_polygon(point, polygon)
-            assert not inside_polygon(point, ex_poly), '%s' %str(point)                        
 
 

inundation/pyvolution/util.py

@@ -5 +5 @@
 """
 
-#FIXME: Move polygon stuff out
+#Soon to be obsoleted here
+
+
+#FIXME: Obsolete this shortcut
+from utilities.numerical_tools import ensure_numeric
+
+
+import utilities.polygon
+
+def point_on_line(*args, **kwargs):
+    """Temporary Interface to new location"""
+
+    print 'point_on_line has moved from util.py.  ',
+    print 'Please use "from utilities.polygon import point_on_line"'
+
+    return utilities.polygon.point_on_line(*args, **kwargs)     
+    
+def inside_polygon(*args, **kwargs):
+    """Temporary Interface to new location"""
+
+    print 'inside_polygon has moved from util.py.  ',
+    print 'Please use "from utilities.polygon import inside_polygon"'
+
+    return utilities.polygon.inside_polygon(*args, **kwargs)    
+    
+def outside_polygon(*args, **kwargs):
+    """Temporary Interface to new location"""
+
+    print 'outside_polygon has moved from util.py.  ',
+    print 'Please use "from utilities.polygon import outside_polygon"'
+
+    return utilities.polygon.outside_polygon(*args, **kwargs)    
+
+
+def separate_points_by_polygon(*args, **kwargs):
+    """Temporary Interface to new location"""
+
+    print 'separate_points_by_polygon has moved from util.py.  ',
+    print 'Please use "from utilities.polygon import separate_points_by_polygon"'
+
+    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)    
+
+
+from utilities.polygon import Polygon_function #No warning
+
+
+def read_polygon(*args, **kwargs):
+    """Temporary Interface to new location"""
+
+    print 'read_polygon has moved from util.py.  ',
+    print 'Please use "from utilities.polygon import read_polygon"'
+
+    return utilities.polygon.read_polygon(*args, **kwargs)    
+
+
+def populate_polygon(*args, **kwargs):
+    """Temporary Interface to new location"""
+
+    print 'populate_polygon has moved from util.py.  ',
+    print 'Please use "from utilities.polygon import populate_polygon"'
+
+    return utilities.polygon.populate_polygon(*args, **kwargs)    
+
+def read_xya(filename):
+    """Read simple xya file, possibly with a header in the first line, just
+    x y [attributes]
+    separated by whitespace
+
+    Format can be either ASCII or NetCDF
+
+    Return list of points, list of attributes and
+    attribute names if present otherwise None
+    """
+    print "read_xya is obsolete.  Use import_points_file from load_mesh.loadASCII"
+    #FIXME: Probably obsoleted by similar function in load_ASCII
+    #FIXME: Name read_data_points (or see 'load' in pylab)
+
+    
+    from load_mesh.loadASCII import import_points_file
+
+    points_dict = import_points_file(filename)
+    return points_dict['pointlist'], points_dict['attributelist']
+
+   
+
+
+#Normal util stuff
 
 def angle(v):
@@ -73 +159 @@
 
 
-def point_on_line(x, y, x0, y0, x1, y1):
-    """Determine whether a point is on a line segment
-
-    Input: x, y, x0, x0, x1, y1: where
-        point is given by x, y
-        line is given by (x0, y0) and (x1, y1)
-
-    """
-
-    from Numeric import array, dot, allclose
-    from math import sqrt
-
-    a = array([x - x0, y - y0])
-    a_normal = array([a[1], -a[0]])
-
-    b = array([x1 - x0, y1 - y0])
-
-    if dot(a_normal, b) == 0:
-        #Point is somewhere on the infinite extension of the line
-
-        len_a = sqrt(sum(a**2))
-        len_b = sqrt(sum(b**2))
-        if dot(a, b) >= 0 and len_a <= len_b:
-           return True
-        else:
-           return False
-    else:
-      return False
-
-def ensure_numeric(A, typecode = None):
-    """Ensure that sequence is a Numeric array.
-    Inputs:
-        A: Sequance. If A is already a Numeric array it will be returned
-                     unaltered
-                     If not, an attempt is made to convert it to a Numeric
-                     array
-        typecode: Numeric type. If specified, use this in the conversion.
-                                If not, let Numeric decide
-
-    This function is necessary as array(A) can cause memory overflow.
-    """
-
-    from Numeric import ArrayType, array
-
-    if typecode is None:
-        if type(A) == ArrayType:
-            return A
-        else:
-            return array(A)
-    else:
-        if type(A) == ArrayType:
-            if A.typecode == typecode:
-                return array(A)  #FIXME: Shouldn't this be just return A?
-            else:
-                return A.astype(typecode)
-        else:
-            return array(A).astype(typecode)
-
-
-
 def file_function(filename,
                   domain = None,
@@ -393 +419 @@
 
 
-def read_xya(filename):
-    """Read simple xya file, possibly with a header in the first line, just
-    x y [attributes]
-    separated by whitespace
-
-    Format can be either ASCII or NetCDF
-
-    Return list of points, list of attributes and
-    attribute names if present otherwise None
-    """
-    print "read_xya is obsolete.  Use import_points_file from load_mesh.loadASCII"
-    #FIXME: Probably obsoleted by similar function in load_ASCII
-    #FIXME: Name read_data_points (or see 'load' in pylab)
-
-    
-    from load_mesh.loadASCII import import_points_file
-
-    points_dict = import_points_file(filename)
-    return points_dict['pointlist'], points_dict['attributelist']
-
-
-#####################################
-#POLYGON STUFF
-#
-#FIXME: All these should be put into new module polygon.py
-
-
-
-#These have been redefined to use separate_by_polygon which will
-#put all inside indices in the first part of the indices array and
-#outside indices in the last
-
-def inside_polygon(points, polygon, closed = True, verbose = False):
-    """See separate_points_by_polygon for documentation
-    """
-
-    from Numeric import array, Float, reshape
-
-    if verbose: print 'Checking input to inside_polygon'
-    try:
-        points = ensure_numeric(points, Float)
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg
-
-    try:
-        polygon = ensure_numeric(polygon, Float)
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg
-
-
-
-    if len(points.shape) == 1:
-        one_point = True
-        points = reshape(points, (1,2))
-    else:
-        one_point = False
-
-    indices, count = separate_points_by_polygon(points, polygon,
-                                                closed, verbose)
-
-    if one_point:
-        return count == 1
-    else:
-        return indices[:count]
-
-def outside_polygon(points, polygon, closed = True, verbose = False):
-    """See separate_points_by_polygon for documentation
-    """
-
-    from Numeric import array, Float, reshape
-
-    if verbose: print 'Checking input to outside_polygon'
-    try:
-        points = ensure_numeric(points, Float)
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg
-
-    try:
-        polygon = ensure_numeric(polygon, Float)
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg
-
-
-
-    if len(points.shape) == 1:
-        one_point = True
-        points = reshape(points, (1,2))
-    else:
-        one_point = False
-
-    indices, count = separate_points_by_polygon(points, polygon,
-                                                closed, verbose)
-
-    if one_point:
-        return count != 1
-    else:
-        return indices[count:][::-1]  #return reversed
-
-
-def separate_points_by_polygon(points, polygon,
-                               closed = True, verbose = False):
-    """Determine whether points are inside or outside a polygon
-
-    Input:
-       points - Tuple of (x, y) coordinates, or list of tuples
-       polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be
-       regarded as belonging to the polygon (closed = True)
-       or not (closed = False)
-
-    Outputs:
-       indices: array of same length as points with indices of points falling
-       inside the polygon listed from the beginning and indices of points
-       falling outside listed from the end.
-
-       count: count of points falling inside the polygon
-
-       The indices of points inside are obtained as indices[:count]
-       The indices of points outside are obtained as indices[count:]
-
-
-    Examples:
-       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
-
-       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
-       will return the indices [0, 2, 1] and count == 2 as only the first
-       and the last point are inside the unit square
-
-    Remarks:
-       The vertices may be listed clockwise or counterclockwise and
-       the first point may optionally be repeated.
-       Polygons do not need to be convex.
-       Polygons can have holes in them and points inside a hole is
-       regarded as being outside the polygon.
-
-    Algorithm is based on work by Darel Finley,
-    http://www.alienryderflex.com/polygon/
-
-    """
-
-    from Numeric import array, Float, reshape, Int, zeros
-
-
-    #Input checks
-    try:
-        points = ensure_numeric(points, Float)
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg
-
-    try:
-        polygon = ensure_numeric(polygon, Float)
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices'
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns'
-
-    assert len(points.shape) == 2,\
-       'Points array must be a 2d array'
-
-    assert points.shape[1] == 2,\
-       'Points array must have two columns'
-
-    N = polygon.shape[0] #Number of vertices in polygon
-    M = points.shape[0]  #Number of points
-
-    px = polygon[:,0]
-    py = polygon[:,1]
-
-    #Used for an optimisation when points are far away from polygon
-    minpx = min(px); maxpx = max(px)
-    minpy = min(py); maxpy = max(py)
-
-
-    #Begin main loop
-    indices = zeros(M, Int)
-
-    inside_index = 0    #Keep track of points inside
-    outside_index = M-1 #Keep track of points outside (starting from end)
-
-    for k in range(M):
-
-        if verbose:
-            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
-
-        #for each point
-        x = points[k, 0]
-        y = points[k, 1]
-
-        inside = False
-
-        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
-            #Check polygon
-            for i in range(N):
-                j = (i+1)%N
-
-                #Check for case where point is contained in line segment
-                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
-                    if closed:
-                        inside = True
-                    else:
-                        inside = False
-                    break
-                else:
-                    #Check if truly inside polygon
-                    if py[i] < y and py[j] >= y or\
-                       py[j] < y and py[i] >= y:
-                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
-                            inside = not inside
-
-        if inside:
-            indices[inside_index] = k
-            inside_index += 1
-        else:
-            indices[outside_index] = k
-            outside_index -= 1
-
-    return indices, inside_index
-
-
-def separate_points_by_polygon_c(points, polygon,
-                                 closed = True, verbose = False):
-    """Determine whether points are inside or outside a polygon
-
-    C-wrapper
-    """
-
-    from Numeric import array, Float, reshape, zeros, Int
-
-
-    if verbose: print 'Checking input to separate_points_by_polygon'
-    #Input checks
-    try:
-        points = ensure_numeric(points, Float)
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg
-
-    #if verbose: print 'Checking input to separate_points_by_polygon 2'
-    try:
-        polygon = ensure_numeric(polygon, Float)
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg
-
-    if verbose: print 'check'
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices'
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns'
-
-    assert len(points.shape) == 2,\
-       'Points array must be a 2d array'
-
-    assert points.shape[1] == 2,\
-       'Points array must have two columns'
-
-    N = polygon.shape[0] #Number of vertices in polygon
-    M = points.shape[0]  #Number of points
-
-    from util_ext import separate_points_by_polygon
-
-    if verbose: print 'Allocating array for indices'
-
-    indices = zeros( M, Int )
-
-    if verbose: print 'Calling C-version of inside poly'
-    count = separate_points_by_polygon(points, polygon, indices,
-                                       int(closed), int(verbose))
-
-    return indices, count
-
-
-
-class Polygon_function:
-    """Create callable object f: x,y -> z, where a,y,z are vectors and
-    where f will return different values depending on whether x,y belongs
-    to specified polygons.
-
-    To instantiate:
-
-       Polygon_function(polygons)
-
-    where polygons is a list of tuples of the form
-
-      [ (P0, v0), (P1, v1), ...]
-
-      with Pi being lists of vertices defining polygons and vi either
-      constants or functions of x,y to be applied to points with the polygon.
-
-    The function takes an optional argument, default which is the value
-    (or function) to used for points not belonging to any polygon.
-    For example:
-
-       Polygon_function(polygons, default = 0.03)
-
-    If omitted the default value will be 0.0
-
-    Note: If two polygons overlap, the one last in the list takes precedence
-
-    FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference???
-
-    """
-
-    def __init__(self, regions, default = 0.0):
-
-        try:
-            len(regions)
-        except:
-            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
-            raise msg
-
-
-        T = regions[0]
-        try:
-            a = len(T)
-        except:
-            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
-            raise msg
-
-        assert a == 2, 'Must have two component each: %s' %T
-
-        self.regions = regions
-        self.default = default
-
-
-    def __call__(self, x, y):
-        from util import inside_polygon
-        from Numeric import ones, Float, concatenate, array, reshape, choose
-
-        x = array(x).astype(Float)
-        y = array(y).astype(Float)
-
-        N = len(x)
-        assert len(y) == N
-
-        points = concatenate( (reshape(x, (N, 1)),
-                               reshape(y, (N, 1))), axis=1 )
-
-        if callable(self.default):
-            z = self.default(x,y)
-        else:
-            z = ones(N, Float) * self.default
-
-        for polygon, value in self.regions:
-            indices = inside_polygon(points, polygon)
-
-            #FIXME: This needs to be vectorised
-            if callable(value):
-                for i in indices:
-                    xx = array([x[i]])
-                    yy = array([y[i]])
-                    z[i] = value(xx, yy)[0]
-            else:
-                for i in indices:
-                    z[i] = value
-
-        return z
-
-def read_polygon(filename):
-    """Read points assumed to form a polygon
-       There must be exactly two numbers in each line
-    """
-
-    #Get polygon
-    fid = open(filename)
-    lines = fid.readlines()
-    fid.close()
-    polygon = []
-    for line in lines:
-        fields = line.split(',')
-        polygon.append( [float(fields[0]), float(fields[1])] )
-
-    return polygon
-
-def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
-    """Populate given polygon with uniformly distributed points.
-
-    Input:
-       polygon - list of vertices of polygon
-       number_of_points - (optional) number of points
-       seed - seed for random number generator (default=None)
-       exclude - list of polygons (inside main polygon) from where points should be excluded 
-
-    Output:
-       points - list of points inside polygon
-
-    Examples:
-       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
-       will return five randomly selected points inside the unit square
-    """
-
-    from random import uniform, seed
-
-    seed(seed)
-
-    points = []
-
-    #Find outer extent of polygon
-    max_x = min_x = polygon[0][0]
-    max_y = min_y = polygon[0][1]
-    for point in polygon[1:]:
-        x = point[0]
-        if x > max_x: max_x = x
-        if x < min_x: min_x = x
-        y = point[1]
-        if y > max_y: max_y = y
-        if y < min_y: min_y = y
-
-
-    while len(points) < number_of_points:
-        x = uniform(min_x, max_x)
-        y = uniform(min_y, max_y)
-
-        append = False
-        if inside_polygon( [x,y], polygon ):
-
-            append = True
-            
-            #Check exclusions
-            if exclude is not None:
-                for ex_poly in exclude:
-                    if inside_polygon( [x,y], ex_poly ):
-                        append = False
-                        
-                    
-        if append is True:                
-            points.append([x,y])                
-
-    return points
-
+
+  
 
 
@@ -914 +500 @@
 import compile
 if compile.can_use_C_extension('util_ext.c'):
-    from util_ext import gradient, gradient2, point_on_line
-    separate_points_by_polygon = separate_points_by_polygon_c
+    from util_ext import gradient, gradient2#, point_on_line
+    #separate_points_by_polygon = separate_points_by_polygon_c
 else:
     gradient = gradient_python

inundation/utilities/polygon.py

@@ -36 +36 @@
 
 
-
-#These have been redefined to use separate_by_polygon which will
-#put all inside indices in the first part of the indices array and
-#outside indices in the last
-
 def inside_polygon(points, polygon, closed = True, verbose = False):
-    """See separate_points_by_polygon for documentation
+    """Determine points inside a polygon
+    
+       Functions inside_polygon and outside_polygon have been defined in
+       terms af separate_by_polygon which will put all inside indices in 
+       the first part of the indices array and outside indices in the last
+        
+       See separate_points_by_polygon for documentation         
     """
 
@@ -77 +78 @@
 
 def outside_polygon(points, polygon, closed = True, verbose = False):
-    """See separate_points_by_polygon for documentation
+    """Determine points outside a polygon
+    
+       Functions inside_polygon and outside_polygon have been defined in
+       terms af separate_by_polygon which will put all inside indices in 
+       the first part of the indices array and outside indices in the last
+        
+       See separate_points_by_polygon for documentation         
     """