Changeset 5189

Timestamp:

Apr 2, 2008, 11:36:59 AM (16 years ago)

Author:

duncan

Message:

adding hacky function, points2polygon

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (6 diffs)
• test_data_manager.py (modified) (3 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -1945 +1945 @@
               %(min(times), max(times), len(times))
         print '  Quantities [SI units]:'
+        # Comment out for reduced memory consumption
         for name in ['stage', 'xmomentum', 'ymomentum']:
             q = fid.variables[name][:].flat
@@ -1960 +1961 @@
     try:
         q = fid.variables[quantity][:] 
+        
+        
     except:
         quantity_dict = {}
@@ -1966 +1969 @@
         #Convert quantity expression to quantities found in sww file    
         q = apply_expression_to_dictionary(quantity, quantity_dict)
-
+    #print "q.shape",q.shape
     if len(q.shape) == 2:
         #q has a time component and needs to be reduced along
@@ -2163 +2166 @@
         fid.close()
         return basename_out
+
 
 #Backwards compatibility
@@ -5849 +5853 @@
     return iterate_over
 
+def points2polygon(points_file,
+                    minimum_triangle_angle=3.0):
+    """
+    WARNING: This function is not fully working.  
+    
+    Function to return a polygon returned from alpha shape, given a points file.
+    
+    WARNING: Alpha shape returns multiple polygons, but this function only returns one polygon.
+    
+    """
+    from anuga.pmesh.mesh import Mesh, importMeshFromFile
+    from anuga.shallow_water import Domain    
+    from anuga.pmesh.mesh_interface import create_mesh_from_regions
+    
+    mesh = importMeshFromFile(points_file)
+    mesh.auto_segment()
+    mesh.exportASCIIsegmentoutlinefile("outline.tsh")
+    mesh2 = importMeshFromFile("outline.tsh")
+    mesh2.generate_mesh(maximum_triangle_area=1000000000, minimum_triangle_angle=minimum_triangle_angle, verbose=False)
+    mesh2.export_mesh_file('outline_meshed.tsh')
+    domain = Domain("outline_meshed.tsh", use_cache = False)
+    polygon =  domain.get_boundary_polygon()
+    return polygon 
 
 #-------------------------------------------------------------
@@ -5858 +5885 @@
     import os 
     os.umask(umask)
-

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -25 +25 @@
 from anuga.coordinate_transforms.geo_reference import Geo_reference, \
      DEFAULT_ZONE
+from anuga.geospatial_data.geospatial_data import Geospatial_data
 
 class Test_Data_Manager(unittest.TestCase):
@@ -2752 +2753 @@
         from Numeric import array, zeros, allclose, Float, concatenate, NewAxis
         from Scientific.IO.NetCDF import NetCDFFile
-        from anuga.geospatial_data.geospatial_data import Geospatial_data
-
         # Used for points that lie outside mesh
         NODATA_value = 1758323
@@ -7391 +7390 @@
         assert (file_line2 == 'goodbye screen catcher\n')
         
+ 
+    def test_points2polygon(self):  
+        att_dict = {}
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[0.0, 0.0]])
+        att_dict['elevation'] = array([10.1, 0.0, 10.4])
+        att_dict['brightness'] = array([10.0, 1.0, 10.4])
+        
+        fileName = tempfile.mktemp(".csv")
+        
+        G = Geospatial_data(pointlist, att_dict)
+        
+        G.export_points_file(fileName)
+        
+        polygon = points2polygon(fileName)
+        
+        # This test may fail if the order changes
+        assert (polygon == [[0.0, 0.0],[1.0, 0.0],[0.0, 1.0]])