Changeset 6261

Timestamp:

Feb 3, 2009, 11:51:59 AM (16 years ago)

Author:

ole

Message:

Looped elevation data conversions

Location:

anuga_work/production/busselton/standardised_version

Files:

• build_busselton.py (modified) (4 diffs)
• project.py (modified) (10 diffs)
• run_busselton.py (modified) (7 diffs)

anuga_work/production/busselton/standardised_version/build_busselton.py

@@ -14 +14 @@
 
 # Standard modules
-from os import sep
-from os.path import dirname, basename
+from os.path import dirname
 from os import mkdir, access, F_OK
-from shutil import copy
-import time
-import sys
 
 # Related major packages
@@ -35 +31 @@
 #------------------------------------------------------------------------------
 
-copy_code_files(project.output_build_time_dir,__file__, 
-               dirname(project.__file__)+sep+ project.__name__+'.py' )
+#copy_code_files(project.output_build_time_dir,__file__, 
+#               dirname(project.__file__)+os.sep+ project.__name__+'.py' )
+#start_screen_catcher(project.output_build_time_dir)
 
-start_screen_catcher(project.output_build_time_dir)
-
-print 'USER:    ', project.user
 
 #-------------------------------------------------------------------------------
@@ -50 +44 @@
 #-------------------------------------------------------------------------------
 
-# FIXME(Ole): Here's the Shark Bay way.
+print 'project.poly_all', project.poly_all
+print 'project.combined_dir_name', project.combined_dir_name
 
-## FIXME (Ole): Clip data by interior regions if possible
+# input elevation directory filenames
+#onshore_in_dir_name = project.onshore_in_dir_name
+#coast_in_dir_name = project.coast_in_dir_name
+#coast_in_dir_name1 = project.coast_in_dir_name1
+#offshore_in_dir_name = project.offshore_in_dir_name
+#offshore_in_dir_name1 = project.offshore_in_dir_name1
+#offshore_in_dir_name2 = project.offshore_in_dir_name2
+#offshore_in_dir_name3 = project.offshore_in_dir_name3
+#offshore_in_dir_name4 = project.offshore_in_dir_name4
+#offshore_in_dir_name5 = project.offshore_in_dir_name5
+
+# output elevation directory filenames
+#onshore_dir_name = project.onshore_dir_name
+#coast_dir_name = project.coast_dir_name
+#coast_dir_name1 = project.coast_dir_name1
+#offshore_dir_name = project.offshore_dir_name
+#offshore_dir_name1 = project.offshore_dir_name1
+#offshore_dir_name2 = project.offshore_dir_name2
+#offshore_dir_name3 = project.offshore_dir_name3
+#offshore_dir_name4 = project.offshore_dir_name4
+#offshore_dir_name5 = project.offshore_dir_name5
+
+
+# Create Geospatial data from ASCII files
+geospatial_data = {}
+for filename in project.ascii_grid_filenames:
+    absolute_filename = project.topographies_in_dir + filename
+    convert_dem_from_ascii2netcdf(absolute_filename,
+                                  basename_out=absolute_filename,
+                                  use_cache=True,
+                                  verbose=True)
+    dem2pts(absolute_filename, use_cache=True, verbose=True)
+
+    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts',
+                                                verbose=True)
+
+# Create Geospatial data from TXT files
+for filename in project.point_filenames:
+    absolute_filename = project.topographies_in_dir + filename
+    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
+                                                verbose=True)
+
+
+#print "creates DEMs from asc data"
+#print onshore_in_dir_name, onshore_dir_name
+#print offshore_in_dir_name5, offshore_dir_name5
 #
-#print 'creating geospatial data objects from asc data (via dem and pts formats)'
-#for filename in project.ascii_grid_filenames:
-#    convert_dem_from_ascii2netcdf(filename,
-#                                  basename_out=filename,
-#                                  use_cache=True, verbose=True)
-#    dem2pts(filename, use_cache=True, verbose=True)#
+#convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
 #
-#    geospatial_data += Geospatial_data(file_name=filename + '.pts',
-#                                       verbose=True)
+#convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)#
 ##
 #
-#print 'creating geospatial data objects from txt data'
-#for filename in project.point_filenames:
-#    geospatial_data += Geospatial_data(file_name=filename + '.txt',
-#                                       verbose=True)
-#
-#
-#print 'clip combined geospatial object by bounding polygon'
-#G = geospatial_data.clip(project.bounding_polygon)
-#
-#
-#print 'export combined geospatial data'
-#if access(project.topographies_dir, F_OK) == 0:
-#    mkdir (project.topographies_dir)
-#G.export_points_file(project.combined_dir_name + '.pts')
-
-
-
-print "project.poly_all", project.poly_all
-print "project.combined_dir_name", project.combined_dir_name
-
-# input elevation directory filenames
-onshore_in_dir_name = project.onshore_in_dir_name
-coast_in_dir_name = project.coast_in_dir_name
-coast_in_dir_name1 = project.coast_in_dir_name1
-offshore_in_dir_name = project.offshore_in_dir_name
-offshore_in_dir_name1 = project.offshore_in_dir_name1
-offshore_in_dir_name2 = project.offshore_in_dir_name2
-offshore_in_dir_name3 = project.offshore_in_dir_name3
-offshore_in_dir_name4 = project.offshore_in_dir_name4
-offshore_in_dir_name5 = project.offshore_in_dir_name5
-
-# output elevation directory filenames
-onshore_dir_name = project.onshore_dir_name
-coast_dir_name = project.coast_dir_name
-coast_dir_name1 = project.coast_dir_name1
-offshore_dir_name = project.offshore_dir_name
-offshore_dir_name1 = project.offshore_dir_name1
-offshore_dir_name2 = project.offshore_dir_name2
-offshore_dir_name3 = project.offshore_dir_name3
-offshore_dir_name4 = project.offshore_dir_name4
-offshore_dir_name5 = project.offshore_dir_name5
-# creates DEM from asc data
-print "creates DEMs from asc data"
-convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)
-
 # creates pts file for onshore DEM
-print "creates pts file for onshore DEM"
-dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
+#print "creates pts file for onshore DEM"
+#dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
+#dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
 
 # create onshore pts files
-print'create Geospatial data1 objects from topographies'
-G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+#print'create Geospatial data1 objects from topographies'
+#G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
 
 # create coastal and offshore pts files
-print'create Geospatial data2 objects from topographies'
-G2 = Geospatial_data(file_name = coast_in_dir_name)
-print'create Geospatial data3 objects from topographies'
-G3 = Geospatial_data(file_name = coast_in_dir_name1)
-print'create Geospatial data4 objects from topographies'
-G_off = Geospatial_data(file_name = offshore_in_dir_name)
-print'create Geospatial data5 objects from topographies'
-G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
-print'create Geospatial data6 objects from topographies'
-G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
-print'create Geospatial data7 objects from topographies'
-G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
-print'create Geospatial data8 objects from topographies'
-G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
-print'create Geospatial data9 objects from topographies'
-G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
+#print'create Geospatial data2 objects from topographies'
+#G2 = Geospatial_data(file_name = coast_in_dir_name)
+#print'create Geospatial data3 objects from topographies'
+#G3 = Geospatial_data(file_name = coast_in_dir_name1)
+#print'create Geospatial data4 objects from topographies'
+#G_off = Geospatial_data(file_name = offshore_in_dir_name)
+#print'create Geospatial data5 objects from topographies'
+#G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
+#print'create Geospatial data6 objects from topographies'
+#G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
+#print'create Geospatial data7 objects from topographies'
+#G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
+#print'create Geospatial data8 objects from topographies'
+#G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
+#print'create Geospatial data9 objects from topographies'
+#G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
 
 
@@ -143 +131 @@
 #-------------------------------------------------------------------------------
 
-print'add all geospatial objects'
-G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
+#print'add all geospatial objects'
+#G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
 
-print'clip combined geospatial object by bounding polygon'
+print 'Add geospatial objects except', project.offshore_name5
+print geospatial_data.keys()
+G = None
+for key in geospatial_data:
+    if key != project.offshore_name5:
+        G += geospatial_data[key]
+        
+print 'Clip (outside) combined geospatial object except %s by area of interest.'
 G_clip = G.clip_outside(project.poly_aoi1)
-G_all = G_clip + G_off5
+
+print 'Clip combined geospatial object by bounding polygon'
+G_all = G_clip + geospatial_data[project.offshore_name5]   # G_off5
 G_clipped = G_all.clip(project.poly_all)
 
-print'export combined DEM file'
-if access(project.topographies_dir,F_OK) == 0:
-    mkdir (project.topographies_dir)
+print 'Export combined DEM file'
+#if access(project.topographies_dir,F_OK) == 0:
+#    mkdir (project.topographies_dir)
 G_clipped.export_points_file(project.combined_dir_name + '.pts')
-G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
+print 'Do txt version too'
+G_clipped.export_points_file(project.combined_dir_name + '.txt') # Use for comparision in ARC
 

anuga_work/production/busselton/standardised_version/project.py

@@ -6 +6 @@
 #------------------------------------------------------------------------------
 
+import os
 from os import sep, environ, getenv, getcwd
 from os.path import expanduser
@@ -52 +53 @@
 finaltime=80000         # final time for simulation
 
-setup='final'  # Final can be replaced with trial or basic.
+setup='trial'  # Final can be replaced with trial or basic.
                # Either will result in a coarser mesh that will allow a
                # faster, but less accurate, simulation.
@@ -58 +59 @@
 if setup =='trial':
     print'trial'
-    res_factor=10
-    time_thinning=48
+    res_factor=100
+    time_thinning=96
     yieldstep=240
 if setup =='basic': 
@@ -102 +103 @@
 #------------------------------------------------------------------------------
 
-# FIXME(Ole): Use lists of elevation data files as in Shark Bay study:
-## elevation data filenames
-#ascii_grid_filenames = ['10m_dem_without_survey', '50m_dem_without_10m_dem',
-#                        'bathysteeppt', 'bathyleft', 'bathyright']
-#point_filenames = ['field_survey_north', 'field_survey_south',
-#                   'clipped_bathymetry_final', 'coast_points_final']
-
-
 
 # elevation data used in build_busselton.py
@@ -140 +133 @@
 landward = 'landward_bounding_polygon.csv'
 
+
 #------------------------------------------------------------------------------
 # Output Elevation Data
@@ -164 +158 @@
 # Location of input and output data
 #------------------------------------------------------------------------------
-# where the input data sits
-onshore_in_dir_name = topographies_in_dir + onshore_name #topo
-coast_in_dir_name = topographies_in_dir + coast_name #coastline
-coast_in_dir_name1 = topographies_in_dir + coast_name1 #beach survey
-offshore_in_dir_name = topographies_in_dir + offshore_name #bathymetry
-offshore_in_dir_name1 = topographies_in_dir + offshore_name1 #bathymetry Charts
-offshore_in_dir_name2 = topographies_in_dir + offshore_name2 #Digitised Fairsheet
-offshore_in_dir_name3 = topographies_in_dir + offshore_name3 #250m
-offshore_in_dir_name4 = topographies_in_dir + offshore_name4 #Bunbury DPI
-offshore_in_dir_name5 = topographies_in_dir + offshore_name5 #Busselton Topo
-
-# where the output data sits
-onshore_dir_name = topographies_dir + onshore_name
-
-coast_dir_name = topographies_dir + coast_name
-coast_dir_name1 = topographies_dir + coast_name1
-
-offshore_dir_name = topographies_dir + offshore_name
-offshore_dir_name1 = topographies_dir + offshore_name1
-offshore_dir_name2 = topographies_dir + offshore_name2
-offshore_dir_name3 = topographies_dir + offshore_name3
-offshore_dir_name4 = topographies_dir + offshore_name4
-offshore_dir_name5 = topographies_dir + offshore_name5
-
-# where the combined elevation file sits
+
+ascii_grid_filenames = [onshore_name,   # Topo
+                        offshore_name5] # Busselton Topo
+
+point_filenames = [coast_name,     # Coastline
+                   coast_name1,    # Beach survey
+                   offshore_name,  # Bathymetry
+                   offshore_name1, # Bathymetry Charts
+                   offshore_name2, # Digitised Fairsheet
+                   offshore_name3, # 250m
+                   offshore_name4] # Bunbury DPI
+
+
+# Where the combined elevation file sits
 combined_dir_name = topographies_dir + combined_name
 combined_smaller_name_dir = topographies_dir + combined_smaller_name
 
-# where the mesh sits (this is created during the run_busselton.py)
+# Where the mesh sits (this is created during the run_busselton.py)
 meshes_dir_name = meshes_dir + scenario_name+'.msh'
 
-# where the boundary ordering files sit (this is used within build_boundary.py)
+# Where the boundary ordering files sit (this is used within build_boundary.py)
 order_filename_dir = boundaries_dir + order_filename
 
-# where the landward points of boundary extent sit (this is used within run_busselton.py)
+# Where the landward points of boundary extent sit (this is used within run_busselton.py)
 landward_dir = boundaries_dir + landward
 
-# where the event sts files sits (this is created during the build_boundary.py)
+# Where the event sts files sits (this is created during the build_boundary.py)
 boundaries_dir_event = boundaries_dir + str(event_number) + sep
 boundaries_dir_mux = muxhome
-
-# where the directory of the output filename sits
+urs_boundary_name = os.path.join(boundaries_dir_event,
+                                 scenario_name)
+
+
+# Where the directory of the output filename sits
 output_build_time_dir = output_dir+build_time+dir_comment+sep   #used for build_busselton.py
 output_run_time_dir = output_dir+run_time+dir_comment+sep       #used for run_busselton.py
 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
 
-#w here the directory of the gauges sit
+# Where the directory of the gauges sit
 gauges_dir_name = gauges_dir + gauge_name       #used for get_timeseries.py
 building_in_dir_name = gauges_dir + building + '.csv'    #used for run_building_inundation.py
@@ -218 +204 @@
 #------------------------------------------------------------------------------
 
-#Land, to set the initial stage/water to be offcoast only
+# Land, to set the initial stage/water to be offcoast only
 poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
 
-#Land, to set the initial stage/water to be offcoast only
+# Land, to set the initial stage/water to be offcoast only
 poly_marina = read_polygon(polygons_dir+'initial_condition_marina.csv')
 
@@ -262 +248 @@
 
     
-trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+trigs_min = number_mesh_triangles(interior_regions, poly_all,
+                                  res_poly_all)
 print 'min estimated number of triangles', trigs_min
    
@@ -270 +257 @@
 #------------------------------------------------------------------------------
 
-# exporting asc grid for Busselton
+# ASCII export grid for Busselton
 xminBusselton = 340000
 xmaxBusselton = 352000
@@ -276 +263 @@
 ymaxBusselton = 6280000
 
-# exporting asc grid for Bunbury
+# ASCII export grid for Bunbury
 xminBunbury = 369000
 xmaxBunbury = 381000

anuga_work/production/busselton/standardised_version/run_busselton.py

@@ -25 +25 @@
 
 # Standard modules
-from os import sep
 import os
-from os.path import dirname, basename
-from os import mkdir, access, F_OK
-from shutil import copy
+from os.path import dirname
 import time
-import sys
 
 # Related major packages
@@ -48 +44 @@
 import project  # Definition of file names and polygons
 
-    
+
 #-----------------------------------------------------------------------
 # Copy scripts to time stamped output directory and capture screen
@@ -56 +52 @@
 #copy script must be before screen_catcher
 #copy_code_files(project.output_run_time_dir, __file__, 
-#         dirname(project.__file__)+sep+ project.__name__+'.py' )
+#         dirname(project.__file__)+os.sep+ project.__name__+'.py' )
 #start_screen_catcher(project.output_run_time_dir, myid, numprocs)
 
@@ -95 +91 @@
 print domain.statistics()
 
-
 domain.set_name(project.scenario_name)
 domain.set_datadir(project.output_run_time_dir) 
@@ -111 +106 @@
                       geo_reference=domain.geo_reference)
 domain.set_quantity('stage', IC, use_cache=True, verbose=True)
-
 domain.set_quantity('friction', project.friction) 
-
 domain.set_quantity('elevation', 
                     filename=project.combined_dir_name+'.pts',
@@ -124 +117 @@
 # Setup boundary conditions 
 #-------------------------------------------------------------------------
-
 print 'Set boundary - available tags:', domain.get_boundary_tags()
-
-boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
 
 Br = Reflective_boundary(domain)
 Bd = Dirichlet_boundary([project.tide,0,0])
-
-
-Bf = Field_boundary(boundary_urs_out+'.sts',
+Bf = Field_boundary(project.urs_boundary_name+'.sts',
                     domain, mean_stage=project.tide,
                     time_thinning=1,
@@ -140 +128 @@
                     use_cache=True,
                     verbose=True)
-
 
 domain.set_boundary({'back': Br,