Changeset 2615

Timestamp:

Mar 28, 2006, 11:00:20 AM (18 years ago)

Author:

nick

Message:

update onslow

Location:

production/onslow_2006

Files:

• project.pyc (modified) (previous)
• project_new.py (moved) (moved from production/onslow_2006/project.py) (4 diffs)
• project_old.py (added)
• run_onslow_new.py (moved) (moved from production/onslow_2006/run_onslow.py) (14 diffs)
• run_onslow_old.py (added)

production/onslow_2006/project_new.py

@@ -9 +9 @@
 
 from pmesh.create_mesh import convert_points_from_latlon_to_utm
+
+from coordinate_transforms.redfearn import degminsec2decimal_degrees
+
+from time import localtime, strftime
+
+from os import getcwd
                 
 #Making assumptions about the location of scenario data
 scenario_dir_name = 'onslow_tsunami_scenario_2006'
 
-
 # 250m data to be provided
 coarsename = 'onsl_bathydem250' # get from Neil/Ingo (DEM or topo data)
-"""
+
 # 30m data to be provided
-finename = 'onshore_30' # get from Neil/Ingo (DEM or topo data) 
-"""
+onshore_name = 'onslow_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data) 
 
-# clipping region for fine elevation data
-west = 240000
-east = 340000
-south = 7570000
-north = 7645000
+offshore_name = 'onslow_offshore_points'
+
+boundary_basename = 'SU-AU'
 
 #swollen/ all data output
 basename = 'source'
+
+codename = 'project.py'
 
 if sys.platform == 'win32':
@@ -36 +40 @@
 
 #Derive subdirectories and filenames
+timedir = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+outputdir = home+sep+scenario_dir_name+sep+'output'+sep+timedir+sep
 meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
 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
 boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
+
+codedir = getcwd()+sep
+                                
+codedirname = codedir + 'project.py'
 
 meshname = meshdir + basename
 
 coarsedemname = datadir + coarsename
-"""
-finedemname = datadir + finename
-combineddemname = datadir + 'onslow_combined_elevation'
-"""
+
+onshore_dem_name = datadir + onshore_name
+
+offshore_dem_name = datadir + offshore_name
+
+combined_dem_name = datadir + 'onslow_combined_elevation'
 
 outputname = outputdir + basename  #Used by post processing
@@ -54 +66 @@
 #!gauge_filename = outputdir + 'onslow_gauges.xya'
 #!gauge_outname = outputdir + 'gauges_max_output.xya'
+
+# clipping region for fine elevation data
+eastingmin = 250000
+eastingmax = 330000
+northingmin = 7580000
+northingmax = 7635000
+
+south = degminsec2decimal_degrees(-22,00,0)
+north = degminsec2decimal_degrees(-21,10,0)
+west = degminsec2decimal_degrees(114,30,0)
+east = degminsec2decimal_degrees(115,30,0)
+
+# region for visualisation
+eminviz = 260000
+emaxviz = 320000
+nminviz = 7590000
+nmaxviz = 7630000
 
 #Georeferencing
@@ -97 +126 @@
 poly_direction = [k0, k1, k2, k3]
 
-#!slump_origin = [385000.0, 6255000.0] #Absolute UTM

production/onslow_2006/run_onslow_new.py

@@ -21 +21 @@
 # Related major packages
 from pyvolution.shallow_water import Domain, Reflective_boundary, \
-                            Dirichlet_boundary, Time_boundary
+                            Dirichlet_boundary, Time_boundary, File_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 geospatial_data import add_points_files
 
 # Application specific imports
@@ -30 +31 @@
 from smf import slump_tsunami  # Function for submarine mudslide 
 
+from shutil import copy
+from os import mkdir, access, F_OK
+
+from geospatial_data import *
 
 #-------------------------------------------------------------------------------
@@ -41 +46 @@
 # filenames
 coarsedemname = project.coarsedemname
-'''
-#finedemname = project.finedemname
-'''
+
+onshore_dem_name = project.onshore_dem_name
+
+offshore_points = project.offshore_dem_name
+
 meshname = project.meshname+'.msh'
 
+source_dir = project.boundarydir
+
+# creates copy of code in output dir
+if access(project.outputdir,F_OK) == 0 :
+    mkdir (project.outputdir)
+copy (project.codedirname, project.outputdir + project.codename)
+copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onslow.py')
+
+
+'''
 # coarse data
 convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
 dem2pts(coarsedemname, use_cache=True, verbose=True)
-'''
+
+
 # fine data (clipping the points file to smaller area)
-convert_dem_from_ascii2netcdf(finedemname, use_cache=True, verbose=True)
-dem2pts(finedemname,
+convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
+dem2pts(onshore_dem_name,
         easting_min=project.eastingmin,
         easting_max=project.eastingmax,
@@ -60 +78 @@
         verbose=True)
 
-
-# combining the coarse and fine data
-combine_rectangular_points_files(project.finedemname + '.pts',
-                                 project.coarsedemname + '.pts',
-                                 project.combineddemname + '.pts')
-'''
-
+        '''
+print'create G1'
+G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
+
+print'create G2'
+G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
+
+print'add G1+G2'
+G = G1 + G2
+
+print'export G'
+G.new_export_points_file(project.combined_dem_name + '.pts')
+
+'''
+add_points_files(
+                  project.offshore_dem_name + '.xya',
+                  project.onshore_dem_name + '.pts',
+                  project.combined_dem_name + '.pts')
+'''
 #-------------------------------------------------------------------------------                                 
 # Create the triangular mesh based on overall clipping polygon with a tagged
@@ -73 +103 @@
 #-------------------------------------------------------------------------------
 
-from pmesh.create_mesh import create_mesh_from_regions
+from pmesh.mesh_interface import create_mesh_from_regions
 
 # original
@@ -80 +110 @@
                     [project.poly_thevenard, interior_res],
                     [project.poly_direction, interior_res]]
-
-#FIXME: Fix caching of this one once the mesh_interface is ready
+                    #[project.testpoly, interior_res]]
+print 'number of interior regions', len(interior_regions)
+
 from caching import cache
 _ = cache(create_mesh_from_regions,
@@ -88 +119 @@
                              'left': [2], 'bottom': [3],
                              'bottomright': [4], 'topright': [5]},
-           'resolution': 100000,
+           'maximum_triangle_area': 100000,
            'filename': meshname,           
-           'interior_regions': interior_regions,
-           'UTM': True,
-           'refzone': project.refzone},
+           'interior_regions': interior_regions},
           verbose = True)
 
@@ -106 +135 @@
 print 'Number of triangles = ', len(domain)
 print 'The extent is ', domain.get_extent()
+print domain.statistics()
 
 domain.set_name(project.basename)
 domain.set_datadir(project.outputdir)
-domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_quantities_to_be_stored(['stage'])
 
 
@@ -136 +166 @@
 domain.set_quantity('stage', tide)
 domain.set_quantity('friction', 0.0) 
+print 'hi and file',project.combined_dem_name + '.pts'
 domain.set_quantity('elevation', 
-#                    filename = project.combineddemname + '.pts',
-                    filename = project.coarsedemname + '.pts',
-                    use_cache = True,
-                    verbose = True
+#                    0.
+#                    filename = project.onshore_dem_name + '.pts',
+                    filename = project.combined_dem_name + '.pts',
+#                    filename = project.offshore_dem_name + '.pts',
+                    use_cache = False,
+                    verbose = True,
+                    alpha = 0.1
                     )
-
+print 'hi1'
 
 #-------------------------------------------------------------------------------                                 
@@ -156 +190 @@
 
 cache(ferret2sww,
-      (source_dir+project.boundary_basename,
-       project.boundary_basename), 
+      (source_dir + project.boundary_basename,
+       source_dir + project.boundary_basename), 
       {'verbose': True,
-       'minlat': south-1,
-       'maxlat': north+1,
-       'minlon': west-1,
-       'maxlon': east+1,
-       'origin': project.mesh_origin,
+# note didn't work with the below
+#       'minlat': south - 1,
+#       'maxlat': north + 1,
+#       'minlon': west - 1,
+#       'maxlon': east + 1,
+       'minlat': south,
+       'maxlat': north,
+       'minlon': west,
+       'maxlon': east,
+#       'origin': project.mesh_origin,
+       'origin': domain.geo_reference.get_origin(),
        'mean_stage': tide,
        'zscale': 1,                 #Enhance tsunami
@@ -174 +214 @@
 print 'Available boundary tags', domain.get_boundary_tags()
 
+Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
+                    domain, verbose = True)
 Br = Reflective_boundary(domain)
 Bd = Dirichlet_boundary([tide,0,0])
@@ -182 +224 @@
                    f=lambda t: [(60<t<480)*6, 0, 0])
 
-domain.set_boundary( {'top': Bw, 'topleft': Br,
+domain.set_boundary( {'top': Bf, 'topleft': Bf,
                              'left': Br, 'bottom': Br,
-                             'bottomright': Br, 'topright': Br} )
+                             'bottomright': Br, 'topright': Bf} )
 
 
@@ -190 +232 @@
 # Evolve system through time
 #-------------------------------------------------------------------------------
-
 import time
 t0 = time.time()
@@ -199 +240 @@
     
 print 'That took %.2f seconds' %(time.time()-t0)
+
+print 'finished'