Changeset 2407

Timestamp:

Feb 14, 2006, 5:19:29 PM (18 years ago)

Author:

ole

Message:

Pushed caching into conversion functions and beautified the Sydney scenario

Files:

• inundation/pyvolution/data_manager.py (modified) (4 diffs)
• inundation/pyvolution/general_mesh.py (modified) (1 diff)
• inundation/pyvolution/pmesh2domain.py (modified) (6 diffs)
• production/sydney_2006/project.py (modified) (9 diffs)
• production/sydney_2006/run_sydney_smf.py (modified) (4 diffs)

inundation/pyvolution/data_manager.py

@@ -1074 +1074 @@
 
 
-def dem2pts(basename_in, basename_out=None, verbose=False,
+def dem2pts(basename_in, basename_out=None, 
             easting_min=None, easting_max=None,
-            northing_min=None, northing_max=None):
+            northing_min=None, northing_max=None,
+            use_cache=False, verbose=False,):
     """Read Digitial Elevation model from the following NetCDF format (.dem)
 
@@ -1093 +1094 @@
     points:  (Nx2) Float array
     elevation: N Float array
+    """
+
+
+
+    kwargs = {'basename_out': basename_out,  
+              'easting_min': easting_min,
+              'easting_max': easting_max,
+              'northing_min': northing_min,
+              'northing_max': northing_max,
+              'verbose': verbose}
+
+    if use_cache is True:
+        from caching import cache
+        result = cache(_dem2pts, basename_in, kwargs,
+                       dependencies = [basename_in + '.dem'],
+                       verbose = verbose)
+
+    else:
+        result = apply(_dem2pts, [basename_in], kwargs)
+        
+    return result
+
+
+def _dem2pts(basename_in, basename_out=None, verbose=False,
+            easting_min=None, easting_max=None,
+            northing_min=None, northing_max=None):
+    """Read Digitial Elevation model from the following NetCDF format (.dem)
+
+    Internal function. See public function dem2pts for details.    
     """
 
@@ -1831 +1861 @@
 
 
-
 def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
-                                  verbose=False):
-    """Read Digitial Elevation model from the following ASCII format (.asc)
+                                  use_cache = False,
+                                  verbose = False):
+    """Read Digitial Elevation model from the following ASCII format (.asc)    
 
     Example:
@@ -1864 +1894 @@
     Yshift        10000000.0000000000
     Parameters
+    """
+
+
+
+    kwargs = {'basename_out': basename_out, 'verbose': verbose}
+
+    if use_cache is True:
+        from caching import cache
+        result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
+                       dependencies = [basename_in + '.asc'],
+                       verbose = verbose)
+
+    else:
+        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
+        
+    return result
+
+    
+    
+
+
+
+def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
+                                  verbose = False):
+    """Read Digitial Elevation model from the following ASCII format (.asc)
+
+    Internal function. See public function convert_dem_from_ascii2netcdf for details.
     """
 

inundation/pyvolution/general_mesh.py

@@ -50 +50 @@
     """
 
+    #FIXME: It would be a good idea to use geospatial data as an alternative
+    #input
     def __init__(self, coordinates, triangles,
                  geo_reference=None):

inundation/pyvolution/pmesh2domain.py

@@ -16 +16 @@
     vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
                         ,tagged_elements_dict, geo_reference = \
-                _pmesh_to_domain(mesh_instance=mesh)
+                        _pmesh_to_domain(mesh_instance=mesh)
 
 
@@ -36 +36 @@
     return domain
 
-def pmesh_to_domain_instance(file_name, DomainClass):
+
+
+def pmesh_to_domain_instance(file_name, DomainClass, use_cache = False, verbose = False):
     """
     Converts a mesh file(.tsh or .msh), to a Domain instance.
@@ -44 +46 @@
     DomainClass is the Class that will be returned.
     It must be a subclass of Domain, with the same interface as domain.
-    
+
+    use_cache: True means that caching is attempted for the computed domain.    
+    """
+
+    if use_cache is True:
+        from caching import cache
+        result = cache(_pmesh_to_domain_instance, (file_name, DomainClass),
+                       dependencies = [file_name],                     
+                       verbose = verbose)
+
+    else:
+        result = apply(_pmesh_to_domain_instance, (file_name, DomainClass))        
+        
+    return result
+
+
+
+
+def _pmesh_to_domain_instance(file_name, DomainClass):
+    """
+    Converts a mesh file(.tsh or .msh), to a Domain instance.
+
+    Internal function. See public interface pmesh_to_domain_instance for details
     """
     
@@ -52 +76 @@
     vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
                         ,tagged_elements_dict, geo_reference = \
-                _pmesh_to_domain(file_name=file_name)
+                        _pmesh_to_domain(file_name=file_name)
 
 
@@ -98 +122 @@
     generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
     generated segment tag list: [S1Tag, S2Tag, ...] (list of ints)
-    triangle list:  [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
+    triangle list:  [(point1, point2, point3),(p5,p4, p1),...] (Tuples of integers)
     triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
     triangle attribute list: [T1att, T2att, ...] (list of strings)
@@ -146 +170 @@
     sides = calc_sides(triangles)
     tag_dict = {}
-    for seg, tag in map(None,mesh_dict['segments'],
-                           mesh_dict['segment_tags']):
+    for seg, tag in map(None, mesh_dict['segments'],
+                        mesh_dict['segment_tags']):
         v1 = seg[0]
         v2 = seg[1]

production/sydney_2006/project.py

@@ -3 +3 @@
 """
 
-from os import sep
+from os import sep, environ
 from os.path import expanduser
+from utilities.polygon import read_polygon
 import sys
+
+from pmesh.create_mesh import convert_points_from_latlon_to_utm
+                
+
+
 
 #Making assumptions about the location of scenario data
@@ -14 +20 @@
 finename = 'bathy_dem25' # get from Neil/Ingo (DEM or topo data) Wed 25 Jan
 
+
+
 # creating easting and northing max and min for fine data - region of interest
 eastingmin = 332090
@@ -29 +37 @@
 
 if sys.platform == 'win32':
-    home = '..\..\..\..\..'     #Sandpit's parent dir
+    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
 else:    
     home = expanduser('~')
+    
 
 
@@ -38 +47 @@
 datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
 outputdir = home+sep+scenario_dir_name+sep+'output'+sep
+polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
 
 meshname = meshdir + basename
@@ -44 +54 @@
 combineddemname = datadir + 'sydneytopo'
 outputname = outputdir + basename  #Used by post processing
+
+#csv file of coastline 50m epsilon belt
+manly_polygonname = polygondir + 'manly_polygon_UTM56_coarse' 
+manly_polygon = read_polygon(manly_polygonname + '.csv')
+#print manly_polygon
+
 gauge_filename = outputdir + 'sydney_gauges.xya'
 gauge_outname = outputdir + 'gauges_max_output.xya'
@@ -51 +67 @@
 
 # define clipping polygon
+south = degminsec2decimal_degrees(-34,05,0)
+north = degminsec2decimal_degrees(-33,33,0)
+west = degminsec2decimal_degrees(151,1,0)
+east = degminsec2decimal_degrees(151,30,0)
+p0 = [south, west]
+p1 = [south, east]
+p2 = [north, east]
+p3 = [north, west]
+    
+polygonall, zone = convert_points_from_latlon_to_utm([p0, p1, p2, p3])
+refzone = zone
+
+print 'Got refzone', refzone
+
 dsouth = degminsec2decimal_degrees(-34,05,0)
 dnorth = degminsec2decimal_degrees(-33,33,0)
@@ -72 +102 @@
 dp8 = [dnorth, dwest]
     
-diffpolygonall = [dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7]
-
+diffpolygonall, zone = convert_points_from_latlon_to_utm([dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7])
+# to put chunk back in
+#diffpolygonall = [dp0, dp1, dp2, dp3, dp4, dp8]
+
+
+
+#Interior regions - the Harbour - take 2
 #Interior regions - the Harbour
 harbour_1x = degminsec2decimal_degrees(-33,51,0)
@@ -122 +157 @@
 k142 = [harbour_15x, harbour_15y]
 
-harbour_polygon_2 = [k02, k112, k122, k12, k22, k62, k72, k82, k102, k42, k52] #worked
+harbour_polygon_2, zone = convert_points_from_latlon_to_utm([k02, k112, k122, k12, k22, k62, k72, k82, k102, k42, k52]) #worked
+assert zone == refzone
+
 
 #Interior region - Botany Bay
+
+#Interior region - Botany Bay - take 2
 bb_1x = degminsec2decimal_degrees(-34,3,0)
 bb_1y = degminsec2decimal_degrees(151,2,30)
@@ -157 +196 @@
 j92 = [bb_10x, bb_10y]
 
-botanybay_polygon_2 = [j92, j12, j22, j62, j82, j72, j42] # worked
-
-# around 42km across from harbour(385000,6255000) 
+botanybay_polygon_2, zone = convert_points_from_latlon_to_utm([j92, j12, j22, j62, j82, j72, j42]) # worked
+
+
+# close to botany bay opening (340000,6236000)
+# x0 = 25964
+# y0 = 11049
+# around 10km from botany bay opening (350000,6236000)
+# x0 = 35964
+# y0 = 11049
+# around 21km from botany bay opening (361000,6236000)
+#x0 = 46964
+#y0 = 11049
+
+# not used for sydney scenario, original interior regions listed though
+# setting up problem area for doing just around the harbour
+hsouth = degminsec2decimal_degrees(-33,54,0)
+hnorth = degminsec2decimal_degrees(-33,48,0)
+hwest = degminsec2decimal_degrees(151,0,0)
+heast = degminsec2decimal_degrees(151,30,0)
+
+hp0 = [hsouth, hwest]
+hp1 = [hsouth, heast]
+hp2 = [hnorth, heast]
+hp3 = [hnorth, hwest]
+polygon_h, zone = convert_points_from_latlon_to_utm([hp0, hp1, hp2, hp3])
+
+#Interior regions - the Harbour - take 1
+harbour_south = degminsec2decimal_degrees(-33,53,0)
+harbour_north = degminsec2decimal_degrees(-33,47,0)
+harbour_west = degminsec2decimal_degrees(151,5,0)
+harbour_east = degminsec2decimal_degrees(151,19,0)
+
+#harbour_south1 = degminsec2decimal_degrees(-33,53,0)
+#harbour_south2 = degminsec2decimal_degrees(-33,52,0)
+#harbour_north1 = degminsec2decimal_degrees(-33,45,0)
+#harbour_north2 = degminsec2decimal_degrees(-33,48,0)
+#harbour_west = degminsec2decimal_degrees(151,5,0)
+#harbour_east = degminsec2decimal_degrees(151,19,0)
+
+k0 = [harbour_south, harbour_west]
+k1 = [harbour_south, harbour_east]
+k2 = [harbour_north, harbour_east]
+k3 = [harbour_north, harbour_west]    
+
+harbour_polygon, zone = convert_points_from_latlon_to_utm([k0, k1, k2, k3])
+
+# setting up problem area for doing just around Botany Bay
+bsouth = degminsec2decimal_degrees(-33,56,0)
+bnorth = degminsec2decimal_degrees(-34,3,0)
+bwest = degminsec2decimal_degrees(151,0,0)
+beast = degminsec2decimal_degrees(151,30,0)
+
+bp0 = [bsouth, bwest]
+bp1 = [bsouth, beast]
+bp2 = [bnorth, beast]
+bp3 = [bnorth, bwest]
+polygon_bb, zone = convert_points_from_latlon_to_utm([bp0, bp1, bp2, bp3])
+
+#Interior region - Botany Bay - take 1
+botanybay_south = degminsec2decimal_degrees(-33,58,0)
+botanybay_north = degminsec2decimal_degrees(-34,1,0)
+botanybay_west = degminsec2decimal_degrees(151,5,0)
+botanybay_east = degminsec2decimal_degrees(151,18,0)
+
+j0 = [botanybay_south, botanybay_west]
+j1 = [botanybay_south, botanybay_east]
+j2 = [botanybay_north, botanybay_east]
+j3 = [botanybay_north, botanybay_west]    
+
+botanybay_polygon, zone = convert_points_from_latlon_to_utm([j0, j1, j2, j3])
+assert zone == refzone
+
 #x0 = 28964 + 42000 
 #y0 = 30049 

production/sydney_2006/run_sydney_smf.py

@@ -6 +6 @@
 
 The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and boundary data obtained from a tsunami simulation done with MOST.
+the elevation data and a simulated submarine landslide.
 
-Ole Nielsen, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
+Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
 """
 
-tide = 0       #Australian Height Datum (mean sea level)
 
+#-------------------------------------------------------------------------------# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
 import os
 import time
 
+# Related major packages
+from pyvolution.shallow_water import Domain, Reflective_boundary
+from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from pyvolution.combine_pts import combine_rectangular_points_files 
+from pyvolution.pmesh2domain import pmesh_to_domain_instance
 
-from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
-     Dirichlet_boundary, Time_boundary, Transmissive_boundary
-from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
-     dem2pts
-from pyvolution.pmesh2domain import pmesh_to_domain_instance
-from pyvolution.combine_pts import combine_rectangular_points_files 
-from caching import cache
-import project
-from math import pi, sin
-from smf import slump_tsunami, slide_tsunami, Double_gaussian
-from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 
+# Application specific imports
+import project                 # Definition of file names and polygons
+from smf import slump_tsunami  # Function for submarine mudslide 
 
-# Data preparation
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
 # Convert ASC 2 DEM 2 PTS using source data and store result in source data
 # Do for coarse and fine data
 # Fine pts file to be clipped to area of interest
+#-------------------------------------------------------------------------------
+
+# filenames
 coarsedemname = project.coarsedemname
 finedemname = project.finedemname
@@ -38 +44 @@
 
 # coarse data
-cache(convert_dem_from_ascii2netcdf, coarsedemname, {'verbose': True},
-      dependencies = [coarsedemname + '.asc'],
-      verbose = True)
-      #evaluate = True)
+convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
+dem2pts(coarsedemname, use_cache=True, verbose=True)
 
-cache(dem2pts, coarsedemname, {'verbose': True},
-      dependencies = [coarsedemname + '.dem'],      
-      verbose = True)
+# fine data (clipping the points file to smaller area)
+convert_dem_from_ascii2netcdf(finedemname, use_cache=True, verbose=True)
+dem2pts(finedemname,
+        easting_min=project.eastingmin,
+        easting_max=project.eastingmax,
+        northing_min=project.northingmin,
+        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
 
-# fine data
-cache(convert_dem_from_ascii2netcdf, finedemname, {'verbose': True},
-      dependencies = [finedemname + '.asc'],
-      verbose = True)
-      #evaluate = True)
-
-# clipping the fine data
-cache(dem2pts, finedemname, {'verbose': True,
-      'easting_min': project.eastingmin,
-      'easting_max': project.eastingmax,
-      'northing_min': project.northingmin,
-      'northing_max': project.northingmax},
-      dependencies = [finedemname + '.dem'],
-      #evaluate = True,
-      verbose = True)
 
 # combining the coarse and fine data
@@ -67 +62 @@
                                  project.coarsedemname + '.pts',
                                  project.combineddemname + '.pts')
-                                 
-# Create Triangular Mesh
-# Overall clipping polygon and interior regions defined in project.py
+
+
+#-------------------------------------------------------------------------------                                 
+# Create the triangular mesh based on overall clipping polygon with a tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------                                 
+
 from pmesh.create_mesh import create_mesh_from_regions
 
-# for whole region
-interior_res = 5000 # maximal area of per triangle
+interior_res = 5000 
 interior_regions = [[project.harbour_polygon_2, interior_res],
-                    [project.botanybay_polygon_2, interior_res]] 
+                    [project.botanybay_polygon_2, interior_res]]
 
-m = cache(create_mesh_from_regions,
-            project.diffpolygonall,
-            {'boundary_tags': {'bottom': [0], 
-                            'right1': [1], 'right0': [2],
-                            'right2': [3], 'top': [4], 'left1': [5],
-                            'left2': [6], 'left3': [7]},
-            'resolution': 100000,
-            'filename': meshname,           
-            'interior_regions': interior_regions},
-            #evaluate=True,    
-            verbose = True)
 
-# Setup domain
-domain = cache(pmesh_to_domain_instance, (meshname, Domain),
-               dependencies = [meshname],                     
-               verbose = True)
+#FIXME: Fix caching of this one once the mesh_interface is ready
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project.diffpolygonall,
+          {'boundary_tags': {'bottom': [0], 
+                             'right1': [1], 'right0': [2],
+                             'right2': [3], 'top': [4], 'left1': [5],
+                             'left2': [6], 'left3': [7]},
+           'resolution': 100000,
+           'filename': meshname,           
+           'interior_regions': interior_regions,
+           'UTM': True,
+           'refzone': project.refzone},
+          verbose = True)
 
-# Bring in elevation data
-domain.set_quantity('elevation',
-                     filename = project.combineddemname + '.pts',
-                     use_cache = True,
-                     verbose = True)
-  
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+
+domain = pmesh_to_domain_instance(meshname, Domain,
+                                  use_cache = True,
+                                  verbose = True)
+
 print 'Number of triangles = ', len(domain)
 print 'The extent is ', domain.get_extent()
@@ -105 +106 @@
 domain.set_name(project.basename)
 domain.set_datadir(project.outputdir)
-domain.store = True
-domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
 
-# Setup Initial conditions
-t = slump_tsunami(length=30000.0,
-                  depth=400.0,
-                  slope=6.0,
-                  thickness=176.0, 
-                  radius=3330,
-                  dphi=0.23,
-                  x0=project.slump_origin[0], 
-                  y0=project.slump_origin[1], 
-                  alpha=0.0, 
-                  domain=domain)
-domain.set_quantity('stage', t)
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+
+tsunami_source = slump_tsunami(length=30000.0,
+                               depth=400.0,
+                               slope=6.0,
+                               thickness=176.0, 
+                               radius=3330,
+                               dphi=0.23,
+                               x0=project.slump_origin[0], 
+                               y0=project.slump_origin[1], 
+                               alpha=0.0, 
+                               domain=domain)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+domain.set_quantity('stage', tsunami_source)
 domain.set_quantity('friction', 0)
-   
-# Setup Boundary Conditions
-print domain.get_boundary_tags()
+domain.set_quantity('elevation',
+                    filename = project.combineddemname + '.pts',
+                    use_cache = True,
+                    verbose = True)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions (all reflective)
+#-------------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
 
 Br = Reflective_boundary(domain)
-Bt = Transmissive_boundary(domain)
-Bd = Dirichlet_boundary([0,0,0])
-# 10 min square wave starting at 1 min, 6m high
-Bw = Time_boundary(domain=domain,
-                   f=lambda t: [(6<t<606)*6, 0, 0])
-
 domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
                       'right2': Br, 'top': Br, 'left1': Br,
                       'left2': Br, 'left3': Br} )
 
-# Evolve
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
 import time
 t0 = time.time()