Changeset 1855 for production/karratha_2005/create_mesh.py

Timestamp:

Sep 28, 2005, 6:00:46 PM (19 years ago)

Author:

ole

Message:

Better interface to meshing for karratha + some pypar work

File:

• production/karratha_2005/create_mesh.py (modified) (4 diffs)

production/karratha_2005/create_mesh.py

@@ -11 +11 @@
 
 
-
-def create_mesh(polygon, tags, resolution, filename = None):
-    """Create mesh from bounding polygon, tags for all segments and resolution.
-
-    Polygon is a list of points in latitude and longitudes - decimal degrees
-
-
-    Tags is a list of symbolic tags
-
-    Resolution is the maximal area
-    
-
-    FIXME - FUTURE:
-      Use class Point_set
-      Take multiple polygons
-
-      Accept deg, min, sec, e.g.
-      [ [(-20,30,55), (116, 20, 00)], ...]
-      
-    
-    """
-
-    from pmesh.mesh import *
-    from pyvolution.coordinate_transforms.redfearn import redfearn
-    
-
-    #Convert to UTM
-    import project
-    refzone = project.refzone  #FIXME
-
-
+def convert_points_from_latlon_to_utm(polygon, refzone):
     points = []
     for point in polygon:
@@ -48 +18 @@
         points.append([easting, northing])
 
-    polygon = points    
+    return points    
 
-    
-    mesh_origin = project.mesh_origin
 
-    geo = Geo_reference(xllcorner = mesh_origin[1],  #From dem
-                        yllcorner = mesh_origin[2],
-                        zone = refzone)
-   
+
+
+
+def create_region(polygon, tags, refzone):
+    """Create dictionary representation of region bounded by given polygon for use with pmesh
+    """
+
 
     #
@@ -77 +48 @@
     #Create tags
     #E.g. ['wall', 'wall', 'ocean', 'wall']
-    segment_tags = [0]*N
-    for key in tags:
-        indices = tags[key]
-        for i in indices:
-            segment_tags[i] = key
+
+    if tags is not None:
+        segment_tags = [0]*N
+        for key in tags:
+            indices = tags[key]
+            for i in indices:
+                segment_tags[i] = key
     
-    dict['segment_tags'] = segment_tags
+        dict['segment_tags'] = segment_tags
 
+
+
+    return dict
+
+
+
+def create_mesh(bounding_polygon, boundary_tags, resolution,
+                filename = None, interior_regions = None):
+    """Create mesh from bounding polygon, tags for all segments and resolution.
+
+    Polygon is a list of points in latitude and longitudes - decimal degrees
+
+
+    Tags is a list of symbolic tags
+
+    Resolution is the maximal area
+
+    Interior_regions is a list of tuples consisting of (polygon, resolution)
+
+    This function does not allow segments to share points - use underlying
+    pmesh functionality for that
+
+    FIXME - FUTURE:
+      Use class Point_set
+      Take multiple polygons
+
+      Accept deg, min, sec, e.g.
+      [ [(-20,30,55), (116, 20, 00)], ...]
+
+
+    FIXME: This whole thing needs to be completely rethought 
+    
+    """
+
+    from pmesh.mesh import *
+    from pyvolution.coordinate_transforms.redfearn import redfearn
+    from pyvolution.util import populate_polygon
+    
+
+    #Make georef
+    #FIXME: Pass in geo or create automatically somehow
+    import project
+    refzone = project.refzone  #FIXME
+    mesh_origin = project.mesh_origin
+    geo = Geo_reference(xllcorner = mesh_origin[1],  #From dem
+                        yllcorner = mesh_origin[2],
+                        zone = refzone)
 
 
@@ -92 +112 @@
     
     m = Mesh(geo_reference=geo)
-        
-    m.addVertsSegs(dict)
 
 
-    #FIXME: Find centroid automatically
-    inner = m.addRegionEN(points[0][0]+1, points[0][1]+1)
+
+    #Convert to UTM
+    bounding_polygon = convert_points_from_latlon_to_utm(bounding_polygon, refzone)    
+    
+    
+    #Do bounding polygon
+    region_dict = create_region(bounding_polygon, boundary_tags, refzone)    
+    m.addVertsSegs(region_dict)
+
+
+    #Find one point inside region automatically
+    if interior_regions is not None:
+        excluded_polygons = []        
+        for P, res in interior_regions:
+            polygon = convert_points_from_latlon_to_utm(P, refzone)            
+            excluded_polygons.append( polygon )
+    else:
+        excluded_polygons = None
+
+    from Numeric import array
+    
+    [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons)
+    print 'I:', inner_point, resolution
+    inner = m.addRegionEN(inner_point[0], inner_point[1])
     inner.setMaxArea(resolution)
+
+
+
+    #Do interior regions
+    if interior_regions is not None:    
+        for P, res in interior_regions:
+            polygon = convert_points_from_latlon_to_utm(P, refzone)                        
+            region_dict = create_region(polygon, None, refzone)
+
+            m.addVertsSegs(region_dict)
+
+            [inner_point] = populate_polygon(polygon, 1)
+            print 'I', inner_point, res
+            X = m.addRegionEN(inner_point[0], inner_point[1])
+            X.setMaxArea(res)
+            
+
     
     m.generateMesh('pzq28.0za1000000a')