Changeset 3072

Timestamp:

Jun 5, 2006, 12:03:41 PM (18 years ago)

Author:

ole

Message:

Made boundary polygon return absolute UTM coordinates as per ticket:81

Location:

inundation/pyvolution

Files:

• general_mesh.py (modified) (10 diffs)
• mesh.py (modified) (8 diffs)
• test_mesh.py (modified) (2 diffs)

inundation/pyvolution/general_mesh.py

@@ +1 @@
+
+from Numeric import concatenate, reshape, take
+from Numeric import array, zeros, Int, Float, sqrt, sum
 
 from coordinate_transforms.geo_reference import Geo_reference
@@ -65 +68 @@
 
         if verbose: print 'General_mesh: Building basic mesh structure' 
-
-        from Numeric import array, zeros, Int, Float, sqrt, sum
 
         self.triangles = array(triangles,Int)
@@ -192 +193 @@
 
 
-    def get_vertex_coordinates(self, obj = False):
+    def get_vertex_coordinates(self, obj=False, absolute=False):
         """Return all vertex coordinates.
         Return all vertex coordinates for all triangles as an Nx6 array
@@ -202 +203 @@
         FIXME - Maybe use something referring to unique vertices?
 
-        
-        """
+        Boolean keyword argument absolute determines whether coordinates
+        are to be made absolute by taking georeference into account        
+        """
+
+        V = self.vertex_coordinates
+        if absolute is True:
+            V0 = self.geo_reference.get_absolute(V[:,0:2])
+            V1 = self.geo_reference.get_absolute(V[:,2:4])
+            V2 = self.geo_reference.get_absolute(V[:,4:6])            
+            
+            V = concatenate( (V0, V1, V2), axis=1 )
+
 
         if obj is True:
-            from Numeric import concatenate, reshape
-            V = self.vertex_coordinates
             #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
 
@@ -213 +222 @@
             return reshape(V, (3*N, 2))
         else:    
-            return self.vertex_coordinates
-
-
-    def get_vertex_coordinate(self, i, j):
+            return V
+
+
+    def get_vertex_coordinate(self, i, j, absolute=False):
         """Return coordinates for vertex j of the i'th triangle.
         Return value is the numeric array slice [x, y]
         """
-        return self.vertex_coordinates[i, 2*j:2*j+2]
+
+        V = self.get_vertex_coordinates(absolute=absolute)
+        return V[i, 2*j:2*j+2]
+    
+        ##return self.vertex_coordinates[i, 2*j:2*j+2]
 
 
@@ -231 +244 @@
         #See quantity.get_vertex_values
         #FIXME (Ole) - oh yes they should
-
-        from Numeric import zeros, Float
 
         N = self.number_of_elements
@@ -251 +262 @@
         """
 
-        from Numeric import take
-
         if (indices ==  None):
             indices = range(len(self))  #len(self)=number of elements
@@ -259 +268 @@
 
     #FIXME - merge these two (get_vertices and get_triangles)
-    def get_triangles(self, obj = False):
+    def get_triangles(self, obj=False):
         """Get connetivity
         Return triangles (triplets of indices into point coordinates)
@@ -269 +278 @@
 
         if obj is True:
-            from Numeric import array, reshape, Int        
             m = len(self)  #Number of triangles
             M = 3*m        #Total number of unique vertices
@@ -320 +328 @@
 
 
-    def get_extent(self):
+    def get_extent(self, absolute=False):
         """Return min and max of all x and y coordinates
-        """
-
-
-        C = self.get_vertex_coordinates()
+
+        Boolean keyword argument absolute determines whether coordinates
+        are to be made absolute by taking georeference into account
+        """
+        
+
+
+        C = self.get_vertex_coordinates(absolute=absolute)
         X = C[:,0:6:2].copy()
         Y = C[:,1:6:2].copy()

inundation/pyvolution/mesh.py

@@ -426 +426 @@
 
 
-    def get_boundary_polygon(self, verbose = False):
+    def get_boundary_polygon(self, verbose=False):
         """Return bounding polygon for mesh (counter clockwise)
 
@@ -432 +432 @@
         If multiple vertex values are present, the algorithm will select the
         path that contains the mesh.
+
+        All points are in absolute UTM coordinates
         """
         
@@ -438 +440 @@
 
 
-        # FIXME (Ole): Make sure all points are absolute UTM    
-        
         # Get mesh extent
-        xmin, xmax, ymin, ymax = self.get_extent() # FIXME: Make Absolute
+        xmin, xmax, ymin, ymax = self.get_extent(absolute=True)
         pmin = ensure_numeric([xmin, ymin])
         pmax = ensure_numeric([xmax, ymax])        
@@ -456 +456 @@
             if edge_id == 2: a = 0; b = 1
 
-            # FIXME: Make Absolute
-            A = self.get_vertex_coordinate(i, a) # Start
-            B = self.get_vertex_coordinate(i, b) # End
+            A = self.get_vertex_coordinate(i, a, absolute=True) # Start
+            B = self.get_vertex_coordinate(i, b, absolute=True) # End
+
 
             # Take the point closest to pmin as starting point
@@ -488 +488 @@
 
 
-
+            
         #Start with smallest point and follow boundary (counter clock wise)
         polygon = [p0]      # Storage for final boundary polygon
-        point_registry = {} # Keep track of storage to avoid multiple runs around boundary
-                            # This will only be the case if there are more than one candidate
-                            # FIXME (Ole): Perhaps we can do away with polygon and use
-                            # only point_registry to save space.
+        point_registry = {} # Keep track of storage to avoid multiple runs around
+                            # boundary. This will only be the case if there are
+                            # more than one candidate.
+                            # FIXME (Ole): Perhaps we can do away with polygon
+                            # and use only point_registry to save space.
 
         point_registry[tuple(p0)] = 0                            
@@ -509 +510 @@
 
                 if verbose:
-                    print 'Point %s has multiple candidates: %s' %(str(p0), candidate_list)
+                    print 'Point %s has multiple candidates: %s'\
+                          %(str(p0), candidate_list)
 
 
@@ -516 +518 @@
                     v_prev = p0 - polygon[-2] # Vector that leads to p0
                 else:
-                    # FIXME (Ole): What do we do if the first point has multiple candidates?
+                    # FIXME (Ole): What do we do if the first point has multiple
+                    # candidates?
                     # Being the lower left corner, perhaps we can use the
                     # vector [1, 0], but I really don't know if this is completely
@@ -530 +533 @@
                     vc = pc-p0  # Candidate vector
                     
-                    # Angle between each candidate and the previous vector in [-pi, pi]
+                    # Angle between each candidate and the previous vector
+                    # in [-pi, pi]
                     ac = angle(vc, v_prev)
                     if ac > pi: ac = pi-ac

inundation/pyvolution/test_mesh.py

@@ -944 +944 @@
                         'FAILED!') 
 
-    def FIXME_test_mesh_get_boundary_polygon_with_georeferencing(self):
+    def test_mesh_get_boundary_polygon_with_georeferencing(self):
      
         # test
@@ -957 +957 @@
 
         relative_points = geo.change_points_geo_ref(absolute_points)
-        #relative_points = absolute_points #uncomment to make test pass!
+
+        #print 'Relative', relative_points
+        #print 'Absolute', absolute_points       
+
         mesh = Mesh(relative_points, vertices, geo_reference=geo)
         boundary_polygon = mesh.get_boundary_polygon()
+
         assert allclose(absolute_points, boundary_polygon)
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon')
+    #suite = unittest.makeSuite(Test_Mesh,'test_mesh_get_boundary_polygon_with_georeferencing')
     suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()