Index: /production/pt_hedland_2006/get_building_inundation.py
===================================================================
--- /production/pt_hedland_2006/get_building_inundation.py	(revision 3270)
+++ /production/pt_hedland_2006/get_building_inundation.py	(revision 3271)
@@ -29,4 +29,5 @@
 
 # Inputs
+#timestampdir = '' # MSL
 #timestampdir = '' # HAT
 timestampdir = '' # LAT
Index: /production/pt_hedland_2006/project.py
===================================================================
--- /production/pt_hedland_2006/project.py	(revision 3270)
+++ /production/pt_hedland_2006/project.py	(revision 3271)
@@ -10,19 +10,20 @@
 from coordinate_transforms.redfearn import degminsec2decimal_degrees
 from time import localtime, strftime
-
 from geospatial_data import *
 
-
-
-
 #Making assumptions about the location of scenario data
+state = 'western_australia'
 scenario_dir_name = 'pt_hedland_tsunami_scenario_2006'
 
 # onshore data from 30m DTED level 2
-onshore_name = 'pt_hedland_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data)
+onshore_name_dted = 'pt_hedland_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data)
+onshore_name_dli = 'pt_hedland_onshore_20m_dli' # get from Neil/Ingo (DEM or topo data)
+
 # offshore data from GA digitised charts
 offshore_name1 = 'pt_hedland_offshore_points'
+
 # offshore data from AHO fairsheets
 offshore_name2 = 'pt_hedland_offshore_points_fairsheet'
+
 # coastline developed from 30m DTED
 coast_name = 'pt_hedland_coastline_points.xya'
@@ -36,11 +37,13 @@
 
 if sys.platform == 'win32':
-    home = getenv('INUNDATIONHOME')
+    home = getenv('INUNDATIONHOME') #Sandpit's parent dir
 #    python_home = getenv('PWD')     
-#    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
     #user = basename(getenv('USERPROFILE'))
     #print 'USER:', user
-else:    
-    home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'1'+sep+'cit'+sep+'risk_assessment_methods_project'+sep+'inundation')     
+else:
+    # original
+    #home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'1'+sep+'cit'+sep+'risk_assessment_methods_project'+sep+'inundation')
+    # update to perlite 2
+    home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'2'+sep+'cit'+sep+'inundation'+sep+'data')     
     user = getenv('LOGNAME')
     print 'USER:', user
@@ -48,21 +51,23 @@
 #Derive subdirectories and filenames
 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
-#print 'home', home
-outputtimedir = home+sep+scenario_dir_name+sep+'output'+sep+time+sep
+outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep+time+sep
+
 #print 'outputtimedir', outputtimedir
-meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
-datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
-gaugedir = home+sep+scenario_dir_name+sep+'gauges'+sep
-polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
-boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
-#output dir without time
-outputdir = home+sep+scenario_dir_name+sep+'output'+sep
-tidedir = home+sep+scenario_dir_name+sep+'tide_data'+sep
+#meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
+#datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
+#gaugedir = home+sep+scenario_dir_name+sep+'gauges'+sep
+#polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
+#boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
+#outputdir = home+sep+scenario_dir_name+sep+'outputs'+sep
+#tidedir = home+sep+scenario_dir_name+sep+'tide_data'+sep
 
-#print'bound', boundarydir
+meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
+datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
+gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+boundarydir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
+tidedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'tide_data'+sep
 
-#gauge_filename = gaugedir + 'onslow_gauges.xya'
-#for MOST
-#gauge_filename = gaugedir + 'pt_hedland_gauges.xya'
 gauge_filename = gaugedir + 'gauge_location_port_hedland.csv'
 buildings_filename = gaugedir + 'pt_hedland_res.csv'
@@ -72,24 +77,9 @@
 coast_filename = datadir + coast_name
 
-# boundary source data
-#MOST_dir = 'f:'+sep+'3'+sep+'ehn'+sep+'users'+sep+'davidb'+sep+'tsunami'+sep+'WA_project'+sep+'SU-AU_90'+sep+'most_2'+sep+'detailed'+sep
-
-#print 'name', __name__
-#print 'path', __file__
-#codedir = getcwd()+sep
-
-#project_code_name = __name__ 
-                           
-#project_code_dir_name = __file__
-
 meshname = meshdir + basename
-
-onshore_dem_name = datadir + onshore_name
-
+onshore_dem_name = datadir + onshore_name_dted
 offshore_dem_name1 = datadir + offshore_name1
 offshore_dem_name2 = datadir + offshore_name2
-
 combined_dem_name = datadir + 'pt_hedland_combined_elevation'
-
 outputname = outputtimedir + basename  #Used by post processing
 
@@ -107,5 +97,4 @@
 
 # region to export (used from export_results.py)
-
 e_min_area = 633000
 e_max_area = 690000
@@ -113,5 +102,5 @@
 n_max_area = 7761000
 
-refzone = 50 # confirm with Hamish
+refzone = 50 
 
 # bounding polygon provided by Hamish
@@ -155,11 +144,2 @@
 
 poly_region = [j0, j1, j2, j3]
-
-coast_buffer_file = datadir+'pts2ascii_test.xya'
-G = Geospatial_data(file_name=coast_buffer_file,delimiter=' ')
-poly_coast = list(G.get_data_points())
-#print 'get_data_points()',G.get_data_points()
-#print 'get_',poly_region
-
-
-
Index: /production/pt_hedland_2006/run_pt_hedland.py
===================================================================
--- /production/pt_hedland_2006/run_pt_hedland.py	(revision 3270)
+++ /production/pt_hedland_2006/run_pt_hedland.py	(revision 3271)
@@ -1,3 +1,3 @@
-"""Script for running a tsunami inundation scenario for Port Hedland, WA, Australia.
+"""Script for running a tsunami inundation scenario for Onslow, WA, Australia.
 
 Source data such as elevation and boundary data is assumed to be available in
@@ -5,21 +5,15 @@
 The output sww file is stored in project.outputtimedir 
 
-The scenario is defined by a triangular mesh created from project.polygon and
-the elevation data.
-
-Ole Nielsen and Duncan Gray, GA - 2005, Nick Bartzis and Jane Sexton, GA - 2006
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
 """
-
-
-#------------------------------------------------------------------------------
-# Import necessary modules
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------# Import necessary modules
+#-------------------------------------------------------------------------------
 
 # Standard modules
-import os 
-
 from os import sep
 from os.path import dirname, basename
-
 import time
 
@@ -30,18 +24,15 @@
 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
-import project                 # Definition of file names and polygons
-from smf import slump_tsunami  # Function for submarine mudslide 
-
 from shutil import copy
 from os import mkdir, access, F_OK
-
 from geospatial_data import *
 import sys
 from pyvolution.util import Screen_Catcher
 
-#------------------------------------------------------------------------------
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+
+#-------------------------------------------------------------------------------
 # Preparation of topographic data
 #
@@ -49,5 +40,5 @@
 # Do for coarse and fine data
 # Fine pts file to be clipped to area of interest
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------
 
 # filenames
@@ -58,10 +49,5 @@
 source_dir = project.boundarydir
 
-
-
-
-#import sys; sys.exit()
-
-# creates copy of code in output dir if dir doesn't exist
+## creates copy of code in output dir if dir doesn't exist
 if access(project.outputtimedir,F_OK) == 0 :
     mkdir (project.outputtimedir)
@@ -69,15 +55,14 @@
 copy (__file__, project.outputtimedir + basename(__file__))
 
-
-#normal screen output is stored in 
+print 'project.outputtimedir',project.outputtimedir
+
+##normal screen output is stored in 
 screen_output_name = project.outputtimedir + "screen_output.txt"
 screen_error_name = project.outputtimedir + "screen_error.txt"
 
 #used to catch screen output to file
-#sys.stdout = Screen_Catcher(screen_output_name)
-#sys.stderr = Screen_Catcher(screen_error_name)
-
-
-'''
+sys.stdout = Screen_Catcher(screen_output_name)
+sys.stderr = Screen_Catcher(screen_error_name)
+
 # fine data (clipping the points file to smaller area)
 # creates DEM from asc data 
@@ -90,38 +75,38 @@
         northing_min=project.northingmin,
         northing_max= project.northingmax,
-        use_cache=True, 
+        use_cache=True,
         verbose=True)
 
-print'create G1'
+print 'create G1'
 G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
 
-print'create G2'
+print 'create G2'
 G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
 
-print'create G3'
+print 'create G3'
 G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
 
-print'add G1+G2+G3'
+print 'add G1+G2+G3'
 G = G1 + G2 + G3
 
-print'export G'
+print 'export G'
 G.export_points_file(project.combined_dem_name + '.pts')
 
-'''
-#------------------------------------------------------------------------------
+
+#-------------------------------------------------------------------------------                                 
 # 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.mesh_interface import create_mesh_from_regions
 
 region_res = 100000
-coast_res = 2500
-pt_hedland_res = 500
+coast_res = 25000
+pt_hedland_res = 5000
 # derive poly_coast from project.coast_name using alpha_shape
+#interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
 interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
-                    [project.poly_coast, coast_res]]
-#                    [project.poly_region, region_res]]
+                    [project.poly_region, region_res]]
 
 print 'number of interior regions', len(interior_regions)
@@ -139,10 +124,10 @@
         count += 1
 
-figname = 'pt_hedland_polys'
-print'pt_hedland_polys'
-#plot_polygons([project.polyAll, project.poly_pt_hedland, project.poly_region],
-plot_polygons([project.poly_coast],
+if sys.platform == 'win32':
+    #figname = project.outputtimedir + 'pt_hedland_polys'
+    figname = 'pt_hedland_polys'
+    plot_polygons([project.polyAll,project.poly_pt_hedland,project.poly_region],
               figname,
-              verbose = True)
+              verbose = True)    
 
 if count == 0:
@@ -152,31 +137,28 @@
     print 'check %s in production directory' %figname
     import sys; sys.exit()
-
+    
 print 'start create mesh from regions'
 from caching import cache
 _ = cache(create_mesh_from_regions,
           project.polyAll,
+#          project.poly_region,
 #          {'boundary_tags': {'right': [0], 'bottomright': [1],
 #                             'bottomleft': [2], 'left': [3], 'top': [4]},
-          {'boundary_tags': {'topright': [0], 'top': [1],'topleft': [2], 'left': [3],
-                             'bottomleft': [4], 'bottomright': [5], 'right': [6]},
-           'maximum_triangle_area': 10000000,
+#          {'boundary_tags': {'topright': [0], 'top': [1],'topleft': [2], 'left': [3],
+#                             'bottomleft': [4], 'bottomright': [5], 'right': [6]},
+          {'boundary_tags': {'topright': [0],'topleft': [1], 'left': [2],
+                             'bottomleft': [3], 'bottomright': [4], 'right': [5]},
+           'maximum_triangle_area': 500000,
            'filename': meshname,           
            'interior_regions': interior_regions},
-          verbose = True)
-
-
-#------------------------------------------------------------------------------
+          verbose = True, evaluate=True)
+
+
+#-------------------------------------------------------------------------------                                 
 # Setup computational domain
-#------------------------------------------------------------------------------
-
-'''
-domain = pmesh_to_domain_instance(meshname, Domain,
-                                  use_cache = False,
-                                  verbose = True)
-'''
-
+#-------------------------------------------------------------------------------                                 
 domain = Domain(meshname, use_cache = False, verbose = True)
 
+print domain.statistics()
 print 'Number of triangles = ', len(domain)
 print 'The extent is ', domain.get_extent()
@@ -187,9 +169,13 @@
 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
 
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------                                 
 # Setup initial conditions
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------
 
 tide = 0.
+#high
+#tide = 3.6
+#low
+#tide = -3.9
 
 domain.set_quantity('stage', tide)
@@ -204,10 +190,7 @@
                     )
 
-#print 'have sent quantities OK - now exiting'
-#import sys; sys.exit()
-
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------                                 
 # Setup boundary conditions (all reflective)
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------
 print 'start ferret2sww'
 # skipped as results in file SU-AU_clipped is correct for all WA
@@ -235,7 +218,6 @@
        'fail_on_NaN': False,
        'inverted_bathymetry': True},
-      #evaluate = True,
-       verbose = True,
-      dependencies= source_dir + project.boundary_basename + '.sww')
+       evaluate = True,
+       verbose = True,)
 
 
@@ -247,42 +229,36 @@
 Bd = Dirichlet_boundary([tide,0,0])
 
-# 7 min square wave starting at 1 min, 6m high
-Bw = Time_boundary(domain = domain,
-                   f=lambda t: [(60<t<480)*6, 0, 0])
-
-#domain.set_boundary( {'top': Bf, 'left': Br, 'bottomleft': Br, 'bottomright': Br, 'right': Br} )
-#
-domain.set_boundary( {'topright': Bf, 'top': Bf,'topleft': Bf, 'left':  Br,
-                             'bottomleft': Br, 'bottomright': Br, 'right': Br})
+domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd,
+                             'bottomleft': Bd, 'bottomright': Bd, 'right': Bd})
                              
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------                                 
 # Evolve system through time
-#------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------
 import time
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 240, finaltime = 7200): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'top')     
-
-for t in domain.evolve(yieldstep = 120, finaltime = 12600
-                       ,skip_initial_step = True): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'top')     
-
-for t in domain.evolve(yieldstep = 60, finaltime = 19800
-                       ,skip_initial_step = True): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'top')     
+for t in domain.evolve(yieldstep = 240, finaltime = 12240): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topleft')     
+
+for t in domain.evolve(yieldstep = 120, finaltime = 15600
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topleft')     
+
+for t in domain.evolve(yieldstep = 60, finaltime = 22020
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topleft')     
     
-for t in domain.evolve(yieldstep = 120, finaltime = 25200
-                       ,skip_initial_step = True): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'top')     
+for t in domain.evolve(yieldstep = 120, finaltime = 27060
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topleft')     
 
 for t in domain.evolve(yieldstep = 240, finaltime = 36000
                        ,skip_initial_step = True): 
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'top')     
+    domain.write_boundary_statistics(tags = 'topleft')     
   
 print 'That took %.2f seconds' %(time.time()-t0)
