Changeset 4147

Timestamp:

Jan 8, 2007, 5:59:47 PM (18 years ago)

Author:

sexton

Message:

(1) updates to Dampier script based on Perth script (2) minor updates to Onslow report

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/util.py (modified) (1 diff)
• anuga_work/production/dampier_2006/build_dampier.py (modified) (4 diffs)
• anuga_work/production/dampier_2006/project.py (modified) (13 diffs)
• anuga_work/production/dampier_2006/run_dampier.py (modified) (12 diffs)
• anuga_work/production/onslow_2006/make_report.py (modified) (1 diff)
• anuga_work/production/onslow_2006/report/modelling_methodology.tex (modified) (1 diff)
• anuga_work/production/onslow_2006/report/onslow_2006_report.tex (modified) (2 diffs)
• anuga_work/production/onslow_2006/report/timeseriesdiscussion.tex (modified) (2 diffs)
• anuga_work/production/onslow_2006/report/tsunami_scenario.tex (modified) (1 diff)
• anuga_work/production/onslow_2006/report_figures/gaugeAirportRunwaystage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBeadonBayeaststage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBeadonBayweststage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBeadonCreekDocksstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBeadonCreekmouthstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBeadonCreeksouthofdockstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBeadonPointLoadingBerthstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeBindiBindiCommunitystage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeCentredamwallstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeDamoverflowstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeHospitalstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeLightTowerstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugePowerStationstage.png (modified) (previous)
• anuga_work/production/onslow_2006/report_figures/gaugeWestofGroynestage.png (modified) (previous)
• anuga_work/production/perth_2006/build_perth.py (modified) (1 diff)
• anuga_work/production/perth_2006/project.py (modified) (2 diffs)
• anuga_work/production/perth_2006/run_perth.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -1030 +1030 @@
         savefig('profilefig')
 
-    depth_axis = axis([time_min/60.0, time_max/60.0, 0, max(max_depths)*1.1])
+    depth_axis = axis([time_min/60.0, time_max/60.0, -0.1, max(max_depths)*1.1])
     stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
-    stage_axis = axis([time_min/60.0, time_max/60.0, -3.0, 3.0])    
     vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
     mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1])  

anuga_work/production/dampier_2006/build_dampier.py

@@ -41 +41 @@
 # output to file
 #------------------------------------------------------------------------------
-
 copy_code_files(project.output_build_time_dir,__file__, 
                dirname(project.__file__)+sep+ project.__name__+'.py' )
@@ -56 +55 @@
 # Fine pts file to be clipped to area of interest
 #-------------------------------------------------------------------------------
+print"project.bounding_polygon",project.bounding_polygon
+print"project.combined_dir_name",project.combined_dir_name
 
-'''
 # topography directory filenames
+onshore_in_dir_name = project.onshore_in_dir_name
+coast_in_dir_name = project.coast_in_dir_name
+island_in_dir_name = project.island_in_dir_name
+island_in_dir_name1 = project.island_in_dir_name1
+island_in_dir_name2 = project.island_in_dir_name2
+island_in_dir_name3 = project.island_in_dir_name3
+offshore_in_dir_name = project.offshore_in_dir_name
+offshore1_in_dir_name = project.offshore1_in_dir_name
+
 onshore_dir_name = project.onshore_dir_name
 coast_dir_name = project.coast_dir_name
-islands_dir_name = project.islands_dir_name
+island_dir_name = project.island_dir_name
+island_dir_name1 = project.island_dir_name1
+island_dir_name2 = project.island_dir_name2
+island_dir_name3 = project.island_dir_name3
 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
-offshore_dir_name6 = project.offshore_dir_name6
-offshore_dir_name7 = project.offshore_dir_name7
-offshore_dir_name8 = project.offshore_dir_name8
-offshore_dir_name9 = project.offshore_dir_name9
-offshore_dir_name10 = project.offshore_dir_name10
-offshore_dir_name11 = project.offshore_dir_name11
-offshore_dir_name12 = project.offshore_dir_name12
-offshore_dir_name13 = project.offshore_dir_name13
-offshore_dir_name14 = project.offshore_dir_name14
 
 # creates DEM from asc data
-convert_dem_from_ascii2netcdf(onshore_dir_name, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(islands_dir_name, use_cache=True, verbose=True)
+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(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(island_in_dir_name1, basename_out=island_dir_name1, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(island_in_dir_name2, basename_out=island_dir_name2, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(island_in_dir_name3, basename_out=island_dir_name3, use_cache=True, verbose=True)
 
 #creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
 dem2pts(onshore_dir_name,
 #        easting_min=project.eastingmin,
@@ -91 +94 @@
         verbose=True)
 
-#creates pts file for islands DEM
-dem2pts(islands_dir_name, use_cache=True, verbose=True)
+#creates pts file for island DEM
+dem2pts(island_dir_name, use_cache=True, verbose=True)
+dem2pts(island_dir_name1, use_cache=True, verbose=True)
+dem2pts(island_dir_name2, use_cache=True, verbose=True)
+dem2pts(island_dir_name3, use_cache=True, verbose=True)
 
-print'create Geospatial data objects from topographies'
-G1 = Geospatial_data(file_name = project.onshore_dir_name + '.pts')
-G2 = Geospatial_data(file_name = project.coast_dir_name + '.xya')
-G3 = Geospatial_data(file_name = project.islands_dir_name + '.pts')
-G_off = Geospatial_data(file_name = project.offshore_dir_name + '.xya')
-G_off1 = Geospatial_data(file_name = project.offshore_dir_name1 + '.xya')
-G_off2 = Geospatial_data(file_name = project.offshore_dir_name2 + '.xya')
-G_off3 = Geospatial_data(file_name = project.offshore_dir_name3 + '.xya')
-G_off4 = Geospatial_data(file_name = project.offshore_dir_name4 + '.xya')
-G_off5 = Geospatial_data(file_name = project.offshore_dir_name5 + '.xya')
-G_off6 = Geospatial_data(file_name = project.offshore_dir_name6 + '.xya')
-G_off7 = Geospatial_data(file_name = project.offshore_dir_name7 + '.xya')
-G_off8 = Geospatial_data(file_name = project.offshore_dir_name8 + '.xya')
-G_off9 = Geospatial_data(file_name = project.offshore_dir_name9 + '.xya')
-G_off10 = Geospatial_data(file_name = project.offshore_dir_name10 + '.xya')
-G_off11 = Geospatial_data(file_name = project.offshore_dir_name11 + '.xya')
-G_off12 = Geospatial_data(file_name = project.offshore_dir_name12 + '.xya')
-G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya')
-G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya')
+print'create Geospatial data1 objects from topographies'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+print'create Geospatial data2 objects from topographies'
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
+print'create Geospatial data3 objects from topographies'
+G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+print'create Geospatial data4 objects from topographies'
+G4 = Geospatial_data(file_name = island_dir_name1 + '.pts')
+print'create Geospatial data5 objects from topographies'
+G5 = Geospatial_data(file_name = island_dir_name2 + '.pts')
+print'create Geospatial data6 objects from topographies'
+G6 = Geospatial_data(file_name = island_dir_name3 + '.pts')
+print'create Geospatial data7 objects from topographies'
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+print'create Geospatial data8 objects from topographies'
+G_off1 = Geospatial_data(file_name = offshore1_in_dir_name + '.xya')
 
+print'add all geospatial objects'
+G = G1 + G2 + G3 + G4 + G5 + G6 + G_off + G_off1 
 
-
-print'clip nw', project.clip_poly_nw
-print'clip e', project.clip_poly_e
-
-print'reading combined_dir_name'
-G_offshore_data = Geospatial_data(file_name = project.topographies_dir+'dampier_combined_elevation_final.pts')
-
-print'reading offshore_dir_name_old'
-G_offshore_old = Geospatial_data(file_name = project.offshore_dir_name_old + '.pts')
-
-
-print 'G_offshore_data_old_nw', 
-G_nw_name = project.topographies_dir + 'nw_old_data'
-G_offshore_nw = G_offshore_old.clip(project.clip_poly_nw)
-G_offshore_nw.export_points_file(G_nw_name + '.pts')
-G_offshore_nw.export_points_file(G_nw_name + '.xya')
-G_nw = array(G_offshore_nw.get_data_points())
-print'shape of arr nw data', G_nw.shape
-print' max and min of array 0',max(G_nw[:,0]),min(G_nw[:,0])
-print' max and min of array 1',max(G_nw[:,1]),min(G_nw[:,1])
-
-print 'G_offshore_data_old_e'#, G_offshore_old_nw.get_data_points()
-G_e_name = project.topographies_dir+'e_old_data'
-G_offshore_e = G_offshore_old.clip(project.clip_poly_e)
-G_offshore_e.export_points_file(G_e_name + '.pts')
-G_offshore_e.export_points_file(G_e_name + '.xya')
-G_e = array(G_offshore_e.get_data_points())
-print'shape of arr e data', G_e.shape
-print' max and min of array 0',max(G_e[:,0]),min(G_e[:,0])
-print' max and min of array 1',max(G_e[:,1]),min(G_e[:,1])
-
-
-print 'G_offshore_data_mid_e'#, G_offshore_old_nw.get_data_points()
-G_mid_e_name = project.topographies_dir+'mid_e_old_data'
-G_offshore_mid_e = G_offshore_old.clip(project.clip_poly_mid_e)
-G_offshore_mid_e.export_points_file(G_mid_e_name + '.pts')
-G_offshore_mid_e.export_points_file(G_mid_e_name + '.xya')
-G_mid_e = array(G_offshore_mid_e.get_data_points())
-print'shape of arr e data', G_mid_e.shape
-print' max and min of array 0',max(G_mid_e[:,0]),min(G_mid_e[:,0])
-print' max and min of array 1',max(G_mid_e[:,1]),min(G_mid_e[:,1])
-
-print 'G_offshore_data_mid_w'#, G_offshore_old_nw.get_data_points()
-G_mid_w_name = project.topographies_dir+'mid_w_old_data'
-G_offshore_mid_w = G_offshore_old.clip(project.clip_poly_mid_w)
-G_offshore_mid_w.export_points_file(G_mid_w_name + '.pts')
-G_offshore_mid_w.export_points_file(G_mid_w_name + '.xya')
-G_mid_w = array(G_offshore_mid_w.get_data_points())
-print'shape of arr e data', G_mid_w.shape
-print' max and min of array 0',max(G_mid_w[:,0]),min(G_mid_w[:,0])
-print' max and min of array 1',max(G_mid_w[:,1]),min(G_mid_w[:,1])
-
-
-print 'G_offshore_data_old_e'#, G_offshore_old_e.get_data_points() 
-print'add all geospatial objects'
-#G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \
-#    + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \
-#    + G_off13 + G_off14
-
-G = G_offshore_data + G_offshore_mid_w + G_offshore_mid_e
-print'shape of arr G data', G.get_data_points().shape
-
-
-#print'clip combined geospatial object by bounding polygon'
+print'clip combined geospatial object by bounding polygon'
 G_clipped = G.clip(project.bounding_polygon)
 #FIXME: add a clip function to pts
-print'shape of clipped data', G_clipped.get_data_points().shape
+#print'shape of clipped data', G_clipped.get_data_points().shape
 
 print'export combined DEM file'
-if access(project.topographies_time_dir,F_OK) == 0:
-    mkdir (project.topographies_time_dir)
-G_clipped.export_points_file(project.combined_time_dir_final_name + '.pts')
-G_clipped.export_points_file(project.combined_time_dir_final_name + '.xya')
-
+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 + '.xya')
 
 '''
@@ -228 +169 @@
        )
 #       dependencies = source_dir + project.boundary_basename + '.sww')
-
-
-
-
-
-
-
-
-

anuga_work/production/dampier_2006/project.py

@@ -2 +2 @@
 """
 
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
 import sys
-from os import sep, environ, getenv, getcwd
-from os.path import expanduser, basename
-
-from anuga.coordinate_transforms.redfearn import\
-     degminsec2decimal_degrees,\
-     convert_from_latlon_to_utm
-
-from time import localtime, strftime, gmtime, ctime
-from anuga.geospatial_data.geospatial_data import *
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
 from anuga.utilities.system_tools import get_user_name
 
@@ -24 +19 @@
 # INUNDATIONHOME is the inundation directory, not the data directory.
 home += sep +'data'
-#----------------------------------
-# Location and naming of scenario data
-#----------------------------------
-state = 'western_australia'
-scenario_name = 'dampier_tsunami'
-scenario_datas_name = 'dampier_tsunami_scenario_2006'  #name of the directory where the data is stored
-#scenario_datas_name = 'karratha_tsunami_scenario_2005' # Tmp location
-
-#mesh_name = 'elevation50m'
-boundaries_name = 'dampier'
-boundaries_source = 'mag_9_corrected'
-#boundaries_source = 'test'
-
-tide = 2.4
-#tide = 0.0
-
-# topography file names
-onshore_name = 'dli_no_islands'
-coast_name = 'DTED_05_Contour'
-islands_name = 'dted_islands'
-offshore_name = 'XY100003902'
-offshore_name1 = 'XY100003903'
-offshore_name2 = 'XY100003951'
-offshore_name3 = 'XY100006321'
-offshore_name4 = 'XY100011756'
-offshore_name5 = 'XY100014243'
-offshore_name6 = 'XY100014244'
-offshore_name7 = 'XY100021081'
-offshore_name8 = 'XY100021082'
-offshore_name9 = 'XY100021083'
-offshore_name10 = 'XY100021085'
-offshore_name11 = 'XY100021086'
-offshore_name12 = 'XY100026309'
-offshore_name13 = 'XY100026338'
-offshore_name14 = 'XYDM83'
-
-offshore_old = 'elevation50m'
-
-combined_name ='dampier_combined_elevation'
-combined_final_name ='dampier_combined_elevation_final'
-
-gauge_name = 'dampier_gauges_up2.csv'
-
-#Derive subdirectories and filenames
-
-meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep
-topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep
-gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
-polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep
-boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
-#outputdir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'output'+sep
-tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep
-
+
+#time stuff
 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
 gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
-#cctime = strftime('%Y%m%d_%H%M%S',ctime()) #gets time for new dir
 build_time = time+'_build'
 run_time = time+'_run'
-
 print 'gtime: ', gtime
-output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
-output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
-output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep
-output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+
+tide = 0.6
+
+#Making assumptions about the location of scenario data
+state = 'western_australia'
+scenario_name = 'dampier'
+scenario = 'dampier_tsunami_scenario_2006'
+
+# onshore data provided by WA DLI
+onshore_name = 'dampier_dli_ext' # original
+#island
+island_name = 'rott_dli_ext' # original
+island_name1 = 'gard_dli_ext'
+island_name2 = 'carnac_island_dted'
+island_name3 = 'penguin_dted'
+
+# AHO + DPI data + colin French coastline
+coast_name = 'waterline'
+offshore_name = 'dampier_bathymetry'
+offshore1_name = 'missing_fairsheets'
+
+#final topo name
+combined_name ='dampier_combined_elevation'
+combined_smaller_name = 'dampier_combined_elevation_smaller'
+
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
 topographies_time_dir = topographies_dir+build_time+sep
-boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep
-boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
-meshes_time_dir = meshes_dir+build_time+sep
-
-#ideas
-#boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep
-
-gauge_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
-gauge_filename = gauge_dir + 'dampier_gauges_up2.csv'
-
-gauges_dir_name = gauges_dir + gauge_name
+
+# input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+island_in_dir_name = topographies_in_dir + island_name
+island_in_dir_name1 = topographies_in_dir + island_name1
+island_in_dir_name2 = topographies_in_dir + island_name2
+island_in_dir_name3 = topographies_in_dir + island_name3
+
+coast_in_dir_name = topographies_in_dir + coast_name
+offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore1_in_dir_name = topographies_in_dir + offshore1_name
 
 onshore_dir_name = topographies_dir + onshore_name
+island_dir_name = topographies_dir + island_name
+island_dir_name1 = topographies_dir + island_name1
+island_dir_name2 = topographies_dir + island_name2
+island_dir_name3 = topographies_dir + island_name3
+
 coast_dir_name = topographies_dir + coast_name
-islands_dir_name = topographies_dir + islands_name
 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
-offshore_dir_name6 = topographies_dir + offshore_name6
-offshore_dir_name7 = topographies_dir + offshore_name7
-offshore_dir_name8 = topographies_dir + offshore_name8
-offshore_dir_name9 = topographies_dir + offshore_name9
-offshore_dir_name10 = topographies_dir + offshore_name10
-offshore_dir_name11 = topographies_dir + offshore_name11
-offshore_dir_name12 = topographies_dir + offshore_name12
-offshore_dir_name13 = topographies_dir + offshore_name13
-offshore_dir_name14 = topographies_dir + offshore_name14
-
-offshore_dir_name_old = topographies_dir + offshore_old
-
-
-
-#output dir
+
+#final topo files
 combined_dir_name = topographies_dir + combined_name
 combined_time_dir_name = topographies_time_dir + combined_name
-combined_time_dir_final_name = topographies_time_dir + combined_final_name
-
-
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
+#combined_time_dir_final_name = topographies_time_dir + combined_final_name
+
+meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
 meshes_dir_name = meshes_dir + scenario_name
-meshes_time_dir_name = meshes_time_dir + scenario_name
-#output_build_time_dir_name = output_build_time_dir + scenario_name  #Used by post processing
+
+polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep
+tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep
+
+boundaries_source = '????'
+#boundaries locations
+boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
+boundaries_in_dir_name = boundaries_in_dir + scenario_name
+boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + scenario_name
+#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
+#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
+
+#output locations
+output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
+output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
+output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
-boundaries_in_dir_name = boundaries_in_dir + boundaries_name
-boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
-boundaries_dir_name = boundaries_dir + boundaries_name
-
-
-
-# Regions
+
+#gauges
+gauge_name = 'dampier.csv'
+gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+
+
+###############################
+# Domain definitions
+###############################
 
 refzone = 50 
 south = degminsec2decimal_degrees(-20,55,0)
 north = degminsec2decimal_degrees(-20,15,0)
-#north = degminsec2decimal_degrees(-19,15,0)
 west = degminsec2decimal_degrees(116,17,0)
 east = degminsec2decimal_degrees(117,10,0)
-
-#only used to clip boundary condition
-#north_boundary = north + 0.02
-#south_boundary = south - 0.02
-#west_boundary = west - 0.02
-#east_boundary = east + 0.02
-
-south_boundary = degminsec2decimal_degrees(-21,0,0)
-#north_boundary = degminsec2decimal_degrees(-19,00,0)
-north_boundary = degminsec2decimal_degrees(-20,10,0)
-#west_boundary = degminsec2decimal_degrees(116,0,0)
-#east_boundary = degminsec2decimal_degrees(118,00,0)
-west_boundary = degminsec2decimal_degrees(116,00,0)
-east_boundary = degminsec2decimal_degrees(117,20,0)
-
 
 p0 = [south, degminsec2decimal_degrees(116,32,0)]
@@ -174 +128 @@
 p8 = [south, east]
 
-bounding_polygon, zone =\
-                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
-#bounding_polygon, zone =\
-#                  convert_from_latlon_to_utm([p1, p2, p3, p4, p5, p6, p7])
-print bounding_polygon
+poly_all, zone = convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
 refzone = zone
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
-print 'poly area', polygon_area(bounding_polygon)/1000000.0
-
-#Interior regions
-
+print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0
+
+res_poly_all = 100000
+
+###############################
+# Interior region definitions
+###############################
+
+poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20_pts.csv')
+res_pos20_neg20 = 20000
+
+poly_dampier = read_polygon(polygons_dir+'dampier_pts.csv')
+res_dampier = 500
+
+poly_karratha = read_polygon(polygons_dir+'karratha_pts.csv')
+res_karratha = 500
+
+poly_delambre = read_polygon(polygons_dir+'delambre_pts.csv')
+res_delambre = 1000
+
+poly_mainisland = read_polygon(polygons_dir+'mainisland_pts.csv')
+res_mainisland = 1000
+
+poly_NWislands = read_polygon(polygons_dir+'NWislands_pts.csv')
+res_NWislands = 1000
+
+plot_polygons([poly_pos20_neg20,poly_dampier,poly_karratha,poly_delambre,polylmainisland,
+               polyNWislands,poly_all],output_run_time_dir + 'poly_pic')
+
+interior_regions = [[poly_pos20_neg20,res_pos20_neg20],[poly_dampier,res_dampier],
+                    [poly_karratha,res_karratha],[poly_delambre,res_delambre],
+                    [poly_mainisland,res_mainisland],[poly_NWislands,res_NWislands]]
+                    
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid - Dampier
+e_min_area = 474000
+e_max_area = 480000
+n_min_area = 7719000
+n_max_area = 7725000
+
+# exporting asc grid - Karratha
+e_min_area = 
+e_max_area = 
+n_min_area = 
+n_max_area = 
+
+"""
+# used in the CIPMA 2006 scenario
 # CIPMA point of interest
 cipma_latitude = -20.588456
 cipma_longitude = 116.771527
-
 
 k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
@@ -198 +197 @@
 assert zone == refzone
 
+poly_facility = read_polygon(polygons_dir+'facility.csv')
+poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')
+poly_interior = read_polygon(polygons_dir+'interior.csv')
+poly_coast = read_polygon(polygons_dir+'coast_final.csv')
+clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')
+clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')
+clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')
+clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
+
+# exporting asc grid
 e_min_area = 474000
 e_max_area = 480000
@@ -203 +212 @@
 n_max_area = 7725000
 
-poly_facility = read_polygon(polygons_dir+'facility.csv')
-
-poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')
-
-poly_interior = read_polygon(polygons_dir+'interior.csv')
-
-poly_coast = read_polygon(polygons_dir+'coast_final.csv')
-
-clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')
-
-clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')
-
-clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')
-
-clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
-
+# used in the original 2005 scenario
 #Interior regions
 karratha_south = degminsec2decimal_degrees(-20,44,0)
@@ -233 +227 @@
 assert zone == refzone
 
-
 #Interior regions
 dampier_south = degminsec2decimal_degrees(-20,40,0)
@@ -248 +241 @@
 assert zone == refzone
 
-
 #Interior regions
 refinery_south = degminsec2decimal_degrees(-20,37,50)
@@ -263 +255 @@
 assert zone == refzone
 
-
 #Interior region around 468899, 7715177:
 #lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
@@ -279 +270 @@
 point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
 assert zone == refzone
-
 
 #Neils areas around interesting points
@@ -297 +287 @@
 assert zone == refzone
 
-
-
-
 neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
                 degminsec2decimal_degrees(116,41,33)]
@@ -314 +301 @@
                                                   neil2_point4])
 assert zone == refzone
-
-
-
-
 
 #Withnell bay
@@ -333 +316 @@
 assert zone == refzone
 
-
-
-
-
 #Larger Withnell bay
 lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
@@ -351 +330 @@
                                                       
 assert zone == refzone
-
-
-
-
+"""

anuga_work/production/dampier_2006/run_dampier.py

@@ -6 +6 @@
 
 The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and a simulated submarine landslide.
+the elevation data and a simulated tsunami generated with URS code.
 
 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
@@ -23 +23 @@
 import sys
 
-
 # Related major packages
 from anuga.shallow_water import Domain
@@ -32 +31 @@
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
-
-from anuga.geospatial_data.geospatial_data import *
 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
@@ -45 +42 @@
 #------------------------------------------------------------------------------
 
+start_screen_catcher(project.output_run_time_dir, myid, numprocs)
 
 # filenames
-
-#build_time = '20061029_231935_build_tide_24'
-#build_time = '20061030_165746_build_tide_24'
-#build_time = '20061102_215532_build_plus_old_data'
-#build_time = '20061103_055258_build'
 build_time = '20061107_063840_build'
-#build_time = '20061025_153643_build_basic'
 bound_time = '20061102_221245_build'
 
 boundaries_name = project.boundaries_name
-#meshes_time_dir_name = project.meshes_time_dir_name+'.msh'
 meshes_dir_name = project.meshes_dir_name+'.msh'
-#source_dir = project.boundarydir
-#boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name
-boundaries_time_dir_name = project.boundaries_time_dir_name
 boundaries_dir_name = project.boundaries_dir_name
+
 tide = project.tide
 
@@ -70 +59 @@
                  dirname(project.__file__)+sep+ project.__name__+'.py' )
 barrier()
-#start_screen_catcher(project.output_run_time_dir, myid, numprocs)
 
 print 'USER: ', project.user
-#sys.exit()
+print 'min triangles', project.trigs_min,
+print 'Note: This is generally about 20% less than the final amount'
+
 #--------------------------------------------------------------------------
 # Create the triangular mesh based on overall clipping polygon with a
@@ -84 +74 @@
     
     print 'start create mesh from regions'
-    interior_regions = [#[project.karratha_polygon, 25000],
-                    [project.poly_coast, 10000],
-                    [project.poly_pipeline, 2000],
-                    [project.poly_facility, 500]]    
-#                    [project.poly_interior, 1000]]       
-    from anuga.utilities.polygon import plot_polygons
-    figname = project.output_run_time_dir + 'poly_pic'
-    plot_polygons([project.poly_coast,project.poly_pipeline,project.poly_facility,
-               project.bounding_polygon],
-               figname,
-               verbose = True)
-#    if access(project.meshes_time_dir,F_OK) == 0:
-#        mkdir(project.meshes_time_dir)
-#    if access(project.meshes_dir,F_OK) == 0:
-#        mkdir(project.meshes_dir)
-    import sys; sys.exit()
-    print 'start create mesh from regions'
-    interior_regions = [#[project.karratha_polygon, 25000],
-#                    [project.cipma_polygon, 1000],
-#                    [project.poly_pipeline, 5000],
-#                    [project.poly_facility, 500]]    
-                    [project.poly_interior, 1000]]    
-#    meshes_dir_name = project.meshes_dir_name + '.msh'
-
-    create_mesh_from_regions(project.bounding_polygon,
-                         boundary_tags={'back': [7, 8], 'side': [0, 6],
-                                        'ocean': [1, 2, 3, 4, 5]},
-                         maximum_triangle_area=100000,
-                         interior_regions=interior_regions,
-#                         filename=meshes_time_dir_name,
-                         filename=meshes_dir_name,
-                         use_cache=True,
-                         verbose=True)
+    create_mesh_from_regions(project.poly_all,
+                             boundary_tags={'back': [7, 8], 'side': [0, 6],
+                                            'ocean': [1, 2, 3, 4, 5]},
+                             maximum_triangle_area=project.res_poly_all,
+                             interior_regions=project.interior_regions,
+                             filename=meshes_dir_name,
+                             use_cache=True,
+                             verbose=True)
 
 # to sync all processors are ready 
@@ -125 +90 @@
 #-------------------------------------------------------------------------
 print 'Setup computational domain'
-#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
 domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
 print domain.statistics()
 
+"""
 print 'starting to create boundary conditions'
 boundaries_in_dir_name = project.boundaries_in_dir_name
 
-from anuga.shallow_water.data_manager import urs2sww, ferret2sww
-
-print 'maxlat=project.south_boundary, minlat=project.north_boundary', project.south_boundary,project.north_boundary
-print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
-print ' maxlon=project.east',project.east
-
-print 'origin: domain.geo_reference.get_origin()',domain.geo_reference.get_origin()
-
-#import sys; sys.exit()
-
-#if access(project.boundaries_time_dir,F_OK) == 0:
-#    mkdir (project.boundaries_time_dir)
+from anuga.shallow_water.data_manager import urs2sw
+
 # put above distribute
 print 'boundary file is: ',boundaries_dir_name
 from caching import cache
 if myid == 0:
-    cache(ferret2sww,
+    cache(urs2sww,
           (boundaries_in_dir_name,
     #       boundaries_time_dir_name), 
@@ -169 +124 @@
            )
 barrier()
-
+"""
 
 #-------------------------------------------------------------------------
@@ -179 +134 @@
 
     domain.set_quantity('stage', tide)
-    domain.set_quantity('friction', 0.0) 
+    domain.set_quantity('friction', 0.01) 
     #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
     print 'Start Set quantity'
 
-
     domain.set_quantity('elevation', 
-                    filename = project.topographies_dir + build_time + sep + project.combined_final_name + '.pts',
+                    filename = project.combined_dir_name + '.txt',
                     use_cache = True,
                     verbose = True,
                     alpha = 0.1)
-    #domain.set_quantity('elevation', -50)
     print 'Finished Set quantity'
 barrier()
@@ -211 +164 @@
 domain.set_store_vertices_uniquely(False)
 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-domain.set_maximum_allowed_speed(0.0) # Allow a little runoff (0.1 is OK)
+domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
 
 #-------------------------------------------------------------------------
@@ -217 +170 @@
 #-------------------------------------------------------------------------
 print 'Available boundary tags', domain.get_boundary_tags()
-
+print 'domain id', id(domain)
 print 'Reading Boundary file'
-#boundariesname = project.boundaries_dir + '20061101_003322_build'+sep+boundaries_name
-#print'boundariesname',boundariesname
-#Bf = File_boundary(boundaries_time_dir_name + '.sww',
-
-#Bf = File_boundary(boundariesname + '.sww',
-
-print 'domain id', id(domain)
-Bf = File_boundary(boundaries_dir_name + '.sww',
-                  domain, time_thinning=5, use_cache=True, verbose=True)
+#Bf = File_boundary(boundaries_dir_name + '.sww',
+#                  domain, time_thinning=5, use_cache=True, verbose=True)
 
 print 'finished reading boundary file'
@@ -246 +192 @@
 t0 = time.time()
 
-#for t in domain.evolve(yieldstep = 60, finaltime = 34000):
+for t in domain.evolve(yieldstep = 60, finaltime = 34000):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'ocean')
+
+#for t in domain.evolve(yieldstep = 120, finaltime = 9000):
 #    domain.write_time()
 #    domain.write_boundary_statistics(tags = 'ocean')
-
-for t in domain.evolve(yieldstep = 120, finaltime = 9000):
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'ocean')
-#    print 'time: ',time.time()
- 
-    if allclose(t, 6000):
-        domain.set_quantity('xmomentum', 0)
-        domain.set_quantity('ymomentum', 0)
-#import sys; sys.exit()
      
-for t in domain.evolve(yieldstep = 60, finaltime = 28800
-                       ,skip_initial_step = True):
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'ocean')    
-
-for t in domain.evolve(yieldstep = 120, finaltime = 34800
-                       ,skip_initial_step = True): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'ocean')    
+#for t in domain.evolve(yieldstep = 60, finaltime = 28800
+#                       ,skip_initial_step = True):
+#    domain.write_time()
+#    domain.write_boundary_statistics(tags = 'ocean')    
+
+#for t in domain.evolve(yieldstep = 120, finaltime = 34800
+#                       ,skip_initial_step = True): 
+#    domain.write_time()
+#    domain.write_boundary_statistics(tags = 'ocean')    
     
 print 'That took %.2f seconds' %(time.time()-t0)

anuga_work/production/onslow_2006/make_report.py

@@ -229 +229 @@
 s = '\\begin{table} \\begin{center} \n'
 fid.write(s)
-s = '\caption{Defined point locations for %s study area.}' %scenario_name
+s = '\caption{Defined point locations for %s study area.}' %scenario_name.title()
 fid.write(s)
 s = """

anuga_work/production/onslow_2006/report/modelling_methodology.tex

@@ -79 +79 @@
 for a range of probabilities (or return periods). As Figure \ref{fig:probonslow}
 shows, for a given probability, a number of events are possible. The resulting
-impact to Onslow would then vary depending on the source of the event.  Further detail
-on the tsunami scenarios are outlined in Section \ref{sec:tsunamiscenario}.
+impact to Onslow would then vary depending on the source of the event.  The
+tsunami scenarios selected for the tsunami risk assessment 
+are discussed in Section \ref{sec:tsunamiscenario}.
 
 % used for the 2005 report when looking at one event

anuga_work/production/onslow_2006/report/onslow_2006_report.tex

@@ -87 +87 @@
  
 \begin{table} \begin{center} 
-\caption{Defined point locations for onslow study area.}
+\caption{Defined point locations for Onslow study area.}
 \label{table:locations}
 \begin{tabular}{|l|l|l|l|}\hline
@@ -158 +158 @@
  
 
+
    \section{Damage modelling inputs}
      \label{sec:damageinputs}

anuga_work/production/onslow_2006/report/timeseriesdiscussion.tex

@@ -7 +7 @@
 case of MSL, this water level will be 0. As the tsunami wave moves
 through this point, the water height may grow and thus the stage will 
-represent the amplitude of the wave. For an onshore location such as the 
-Light Tower, the actual water depth will be the difference between
-the stage and the elevation at that point. Therefore, at the beginning
-of the simulation, there will be no water onshore and therefore
-the stage and the elevation will be identical.}. Both stage and speed 
+represent the amplitude of the wave.} For an onshore location such as the 
+Light Tower, the actual water depth will be shown rather than the stage. 
+Both stage and speed 
 (in metres/second) for
 each scenario (HAT, MSL and LAT) are shown
@@ -62 +60 @@
 the actual amplitude is the difference between the stage value
 and the initial water level; 2.3 - 1.5}.
-The drawdown of around 4.3 m (i.e. 2.3 - -2) then occurs at around 230 mins
+The drawdown of 4.05 m (i.e. 2.3 - -1.75) then occurs at around 230 mins
 (i.e. 3.8 hours after the event has been generated), before
 the second wave arrives 

anuga_work/production/onslow_2006/report/tsunami_scenario.tex

@@ -1 @@
-The tsunamigenic event used in this report was developed for a 
-preliminary tsunami hazard assessment study delivered by GA
-to FESA in September 2005
-\cite{BC:FESA}. In the assessment, a suite of Mw 9 earthquakes
-were evenly spaced along the Sunda Arc subduction zone and there
-was no consideration of the likelihood of each event. 
-Other less likely sources were not considered, such
-as intra-plate earthquakes near the WA coast, volcanoes, landslides
-or asteroids. 
-In the preliminary assessment,
-the maximum magnitude of earthquakes off Java was considered to be
-at least 8.5 and could potentially be as high as 9.
+% for original scenario
+%The tsunamigenic event used in this report was developed for a 
+%preliminary tsunami hazard assessment study delivered by GA
+%to FESA in September 2005
+%\cite{BC:FESA}. In the assessment, a suite of Mw 9 earthquakes
+%were evenly spaced along the Sunda Arc subduction zone and there
+%was no consideration of the likelihood of each event. 
+%Other less likely sources were not considered, such
+%as intra-plate earthquakes near the WA coast, volcanoes, landslides
+%or asteroids. 
+%In the preliminary assessment,
+%the maximum magnitude of earthquakes off Java was considered to be
+%at least 8.5 and could potentially be as high as 9.
 
-FESA is interested in the ``most frequent worst case scenario''. Whilst
-we currently cannot determine exactly what that event may be, the Mw 9 event 
-provides a plausible worst case scenario. To understand the 
-frequency of these tsunami-genic events,
-GA is building probabilistic
-models to develop a more complete tsunami hazard assessment
-for the Sunda Arc subduction zone,
-due for completion in late 2006. In the preliminary assessment for
-example, it was suggested that while Mw 7 and 8 earthquakes are expected
-to occur with a greater frequency than Mw 9 events, 
-they are likely to pose a comparatively low and more localised hazard to WA. 
+%FESA is interested in the ``most frequent worst case scenario''. Whilst
+%we currently cannot determine exactly what that event may be, the Mw 9 event 
+%provides a plausible worst case scenario. To understand the 
+%frequency of these tsunami-genic events,
+%GA is building probabilistic
+%models to develop a more complete tsunami hazard assessment
+%for the Sunda Arc subduction zone,
+%due for completion in late 2006. In the preliminary assessment for
+%example, it was suggested that while Mw 7 and 8 earthquakes are expected
+%to occur with a greater frequency than Mw 9 events, 
+%they are likely to pose a comparatively low and more localised hazard to WA. 
 
-Figure \ref{fig:mw9} shows the maximum wave height of a tsunami initiated 
-by a Mw 9 event off
-the coast of Java. This event provides the source and
-boundary condition to the
-inundation model presented in Section \ref{sec:anuga}. 
+%Figure \ref{fig:mw9} shows the maximum wave height of a tsunami initiated 
+%by a Mw 9 event off
+%the coast of Java. This event provides the source and
+%boundary condition to the
+%inundation model presented in Section \ref{sec:anuga}. 
 
 
-\begin{figure}[hbt]
+%\begin{figure}[hbt]
 
-  \centerline{ \includegraphics[width=140mm, height=100mm]
-{../report_figures/mw9.jpg}}
+%  \centerline{ \includegraphics[width=140mm, height=100mm]
+%{../report_figures/mw9.jpg}}
 
-  \caption{Maximum wave height (in cms) for a Mw 9 event off the 
-coast of Java}
-  \label{fig:mw9}
-\end{figure}
+%  \caption{Maximum wave height (in cms) for a Mw 9 event off the 
+%coast of Java}
+%  \label{fig:mw9}
+%\end{figure}

anuga_work/production/perth_2006/build_perth.py

@@ -1 @@
-"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+"""Script for running tsunami inundation scenario for Perth, WA, Australia.
 
 Source data such as elevation and boundary data is assumed to be available in

anuga_work/production/perth_2006/project.py

@@ -43 +43 @@
 island_name3 = 'penguin_dted'
 
-# AHO + DPI data
+# AHO + DPI data + colin French coastline
 coast_name = 'waterline'
 offshore_name = 'perth_bathymetry'
@@ -57 +57 @@
 topographies_time_dir = topographies_dir+build_time+sep
 
-#input topo file location
+# input topo file location
 onshore_in_dir_name = topographies_in_dir + onshore_name
 island_in_dir_name = topographies_in_dir + island_name

anuga_work/production/perth_2006/run_perth.py

@@ -1 @@
-"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+"""Script for running tsunami inundation scenario for Perth, WA, Australia.
 
 Source data such as elevation and boundary data is assumed to be available in
@@ -6 +6 @@
 
 The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and a simulated submarine landslide.
+the elevation data and a simulated tsunami generated with URS code.
 
 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006