Changeset 664 for inundation/ga/storm_surge/pyvolution/util.py

Timestamp:

Dec 2, 2004, 3:28:33 PM (20 years ago)

Author:

ole

Message:

Added arbitrary functions to Polygon_function

File:

• inundation/ga/storm_surge/pyvolution/util.py (modified) (7 diffs)

inundation/ga/storm_surge/pyvolution/util.py

@@ -57 +57 @@
     """
     
-    #FIXME: Could do with some optimisation         
+    #FIXME: Could do with some C-optimisation         
         
     from Numeric import array, dot, allclose
@@ -94 +94 @@
        point - 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) 
     
     Output:
@@ -109 +112 @@
        
     Remarks:
-       <Copy from source Franklins code >    
+       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.          
                
-               
-       points at a vertex are considered to be inside          
        
-    Algorithm due to Darel Finley, http://www.alienryderflex.com/polygon/ 
-    
-    
-    FIXME: More comments about algoritm
-    FIXME: Deal with case where point on border of polygon (closed)
-       
+    Algorithm is based on work by Darel Finley, 
+    http://www.alienryderflex.com/polygon/ 
+    
     """    
 
@@ -166 +168 @@
         y = points[k, 1]        
 
-        
         inside = False
         for i in range(N):
@@ -187 +188 @@
         if inside: indices.append(k)
 
-        #print x, y, polygon, inside
-
-    
     if one_point: 
         return inside
@@ -222 +220 @@
     Note: If two polygons overlap, the one last in the list takes precedence
 
-    Note: default can only be a constant.    
     """
     
@@ -253 +250 @@
         x = array(x).astype(Float)
         y = array(y).astype(Float)      
-        x = reshape(x, (len(x), 1))
-        y = reshape(y, (len(y), 1))     
-                
-        points = concatenate( (x,y), axis=1 )
-        
-        z = ones(x.shape[0], Float) * self.default
+        
+        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)
-            for i in indices:
-                z[i] = value 
-
-        #from sys import exit; exit()   
+            
+            #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