Index: /anuga_work/production/MOST_example/cripple_pts_file.py
===================================================================
--- /anuga_work/production/MOST_example/cripple_pts_file.py	(revision 3606)
+++ /anuga_work/production/MOST_example/cripple_pts_file.py	(revision 3606)
@@ -0,0 +1,5 @@
+from anuga.load_mesh.loadASCII import reduce_pts, produce_half_point_files
+
+produce_half_point_files('pt_hedland_combined_elevation.pts',
+                         1000,'_', verbose=True)
+#import load_mesh.loadASCII
Index: /anuga_work/production/MOST_example/project.py
===================================================================
--- /anuga_work/production/MOST_example/project.py	(revision 3606)
+++ /anuga_work/production/MOST_example/project.py	(revision 3606)
@@ -0,0 +1,181 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+Also includes origin for slump scenario.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser, basename
+#from anuga.utilities.polygon import read_polygon
+import sys
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
+from time import localtime, strftime
+from anuga.geospatial_data.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_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_orig'
+
+# offshore data from AHO fairsheets
+offshore_name2 = 'pt_hedland_offshore_points_fairsheet_orig'
+
+# coastline developed from aerial photography and 1.5m DLI contour
+#coast_name = 'pt_hedland_coastline_points_dted2_new'
+coast_name = 'coast_with_extracted_z'
+
+
+boundary_basename = 'SU-AU_clip'
+
+#swollen/ all data output
+basename = 'source'
+
+codename = 'project.py'
+
+if sys.platform == 'win32':
+    home = getenv('INUNDATIONHOME') #Sandpit's parent dir
+#    python_home = getenv('PWD')     
+    user = getenv('USERPROFILE')
+    #print 'USER:', user
+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')
+    
+#Derive subdirectories and filenames
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+outputtimedir = 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
+#outputdir = home+sep+scenario_dir_name+sep+'outputs'+sep
+#tidedir = home+sep+scenario_dir_name+sep+'tide_data'+sep
+
+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 + 'gauge_location_port_hedland.csv'
+gauge_checking = gaugedir + 'gauge_checking_test.csv'
+buildings_filename = gaugedir + 'pt_hedland_res.csv'
+buildings_filename_out = 'pt_hedland_res_modified.csv'
+buildings_filename_damage_out = 'pt_hedland_res_modified_damage.csv'
+community_filename = gaugedir + 'CHINS_v2.csv'
+community_scenario = gaugedir + 'community_pt_hedland.csv'
+tidal_filename = tidedir + 'pt_hedland_tide.txt'
+
+meshname = basename+time
+mesh_elevname = basename+time+'_elev'
+#onshore_dem_name = datadir + onshore_name_dted
+onshore_dem_name = datadir + onshore_name_dli
+offshore_dem_name1 = datadir + offshore_name1
+offshore_dem_name2 = datadir + offshore_name2
+coast_dem_name = datadir + coast_name
+combined_dem_name = 'pt_hedland_combined_elevation'
+outputname = outputtimedir + basename  #Used by post processing
+
+# clipping region to make DEM (pts file) from onshore data
+eastingmin = 594000
+eastingmax = 715000
+northingmin = 7720000
+northingmax = 7880000
+
+# for ferret2sww
+south = degminsec2decimal_degrees(-20,30,0)
+north = degminsec2decimal_degrees(-17,10,0)
+west = degminsec2decimal_degrees(117,00,0)
+east = degminsec2decimal_degrees(120,00,0)
+
+# region to export (used from export_results.py)
+e_min_area = 648000
+e_max_area = 675000
+n_min_area = 7745000
+n_max_area = 7761000
+
+export_region = [[e_min_area,n_min_area],[e_min_area,n_max_area],[e_max_area,n_max_area],[e_max_area,n_min_area]]
+                 
+refzone = 50 
+
+from anuga.coordinate_transforms.redfearn import redfearn
+# boundary up to 50 m contour
+lat1_50 = degminsec2decimal_degrees(-19,20,0)
+lat2_50 = degminsec2decimal_degrees(-19,30,0)
+lat3_50 = degminsec2decimal_degrees(-19,45,0)
+lon1_50 = degminsec2decimal_degrees(119,05,0)
+lon2_50 = degminsec2decimal_degrees(118,20,0)
+lon3_50 = degminsec2decimal_degrees(117,45,0)
+z, easting, northing = redfearn(lat1_50, lon1_50)
+d0_50 = [easting, northing]
+z, easting, northing = redfearn(lat2_50, lon2_50)
+d1_50 = [easting, northing]
+z, easting, northing= redfearn(lat3_50, lon3_50)
+d2_50 = [easting, northing]
+
+d4_50 = [285000, 7585000]
+d6_50 = [330000, 7605000]
+#bounding_poly50 = [p0_50, p1_50, p2_50, d6_50, d5, d4_50]
+
+d0 = [763852.0, 7934358.0]
+d1 = [710987.0, 7925797.0]
+d2 = [658264.0, 7926314.0]
+d3 = [552686.0, 7871580.0]
+#d4 = [604415.81, 7733013.56]
+d4 = [638000.0, 7733013.56]
+#d5 = [656561.15, 7732615.11]
+d5 = [662000.0, 7732615.11]
+#d6 = [708940.32, 7750510.33]
+d6 = [690000.0, 7740510.33]
+#polyAll = [d0, d1, d2, d3, d4, d5, d6]
+#polyAll = [d0_50, d1_50, d2_50, d4, d5, d6]
+# from Hamish
+h0=[629262.17, 7747205.47]
+h1=[552686.00, 7871579.99] #d3
+h2=[658264.00, 7926314.00] #d2
+h3=[710986.99, 7925796.99] #d1
+h4=[763851.99, 7934357.99] #d0
+h5=[701485.21, 7770656.86]
+h6=[698273.75, 7762227.38]
+h7=[698194.23, 7762018.65]
+h8=[691627.41, 7744781.98]
+h9=[679220.75, 7743604.59]
+h10=[653512.59, 7740528.56]
+h11=[634777.71, 7738247.17]
+h12=[629443.86, 7746910.37]
+h13=[629396.84, 7746986.75]
+h14=[629352.32, 7747059.06]
+h15=[629276.24, 7747182.63]
+h16=[629262.17, 7747205.47] #repeat of h0
+# using Hamish's new bounding polygon
+#polyAll = [d0_50, d1_50, d2_50, h16,h15,h14,h13,h12,h11,h10,h9,h8,h7,h6,h5]
+polyAll = [d0_50, d1_50, d2_50, h16,h11,h8,h6, h5]
+
+#Interior region - Pt Hedland town
+i0 = [668000, 7757000]
+i1 = [659000, 7755000]
+i2 = [660000, 7749000]
+i3 = [667000, 7746000]
+i4 = [678000, 7751000]
+
+poly_pt_hedland = [i0, i1, i2, i3, i4]
+
+#Are there other significant features?
+j0 = [670000, 7760000]
+j1 = [633000, 7745000]
+j2 = [665000, 7743000]
+j3 = [690000, 7755000]
+
+poly_region = [j0, j1, j2, j3]
Index: /anuga_work/production/MOST_example/run_pt_hedland.py
===================================================================
--- /anuga_work/production/MOST_example/run_pt_hedland.py	(revision 3606)
+++ /anuga_work/production/MOST_example/run_pt_hedland.py	(revision 3606)
@@ -0,0 +1,215 @@
+"""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
+directories specified by project.py
+The output sww file is stored in project.outputtimedir 
+
+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
+#-------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+import time
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary, Time_boundary, File_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \
+     dem2pts, ferret2sww
+from shutil import copy
+from os import mkdir, access, F_OK
+from anuga.geospatial_data.geospatial_data import *
+import sys
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
+
+from anuga.fit_interpolate.fit import fit_to_mesh_file
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+
+# creates copy of code in output dir if dir doesn't exist
+if access(project.outputtimedir,F_OK) == 0 :
+    mkdir (project.outputtimedir)
+copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
+copy (__file__, project.outputtimedir + basename(__file__))
+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)
+print 'USER:    ', project.user
+
+#-------------------------------------------------------------------------------
+# 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
+meshname = project.meshname + '.tsh'
+mesh_elevname = project.mesh_elevname  + '.tsh'
+source_dir = project.boundarydir
+
+
+#-------------------------------------------------------------------------------                                 
+# 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 anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+region_res = 5000000
+coast_res = 500000
+pt_hedland_res = 500000
+interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
+                    [project.poly_region, region_res]]
+
+print 'number of interior regions', len(interior_regions)
+
+print 'start create mesh from regions'
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project.polyAll,
+          {'boundary_tags': {'topright': [0], 'topleft': [1],
+                             'left': [2], 'bottom0': [3],
+                             'bottom1': [4], 'bottom2': [5],
+			     'bottom3': [6], 'right': [7]},
+           'maximum_triangle_area': 5000000,
+           'filename': meshname,           
+           'interior_regions': interior_regions},
+          verbose = True, evaluate=True)
+
+cache(fit_to_mesh_file,(meshname,
+                 'pt_hedland_combined_elevation_31204' + '.pts',
+                 mesh_elevname),
+      {'verbose': True}
+      #,evaluate = True      
+      ,verbose = False
+      )
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+domain = Domain(mesh_elevname, use_cache = False, verbose = True)
+
+print domain.statistics()
+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.outputtimedir)
+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)
+domain.set_quantity('friction', 0.0) 
+print 'hi and file',project.combined_dem_name + '.pts'
+
+#domain.set_quantity('elevation', 
+#                    filename = project.combined_dem_name + '.pts',
+#                    use_cache = True,
+#                    verbose = True,
+#                    alpha = 0.1
+ #                   )
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions (all reflective)
+#-------------------------------------------------------------------------------
+print 'start ferret2sww'
+# skipped as results in file SU-AU_clipped is correct for all WA
+
+
+south = project.south 
+north = project.north
+west = project.west
+east = project.east
+
+#note only need to do when an SWW file for the MOST boundary doesn't exist
+cache(ferret2sww,
+      (project.boundary_basename,
+       project.boundary_basename+'_'+project.basename), 
+      {'verbose': True,
+       'minlat': south,
+       'maxlat': north,
+       'minlon': west,
+       'maxlon': east,
+#       'origin': project.mesh_origin,
+       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': tide,
+       'zscale': 10,                 #Enhance tsunami
+       'fail_on_NaN': False,
+       'inverted_bathymetry': True},
+       evaluate = True,
+       verbose = True,
+      dependencies = source_dir + project.boundary_basename + '.sww')
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Bf = File_boundary(project.boundary_basename+'_'+project.basename + '.sww', 
+                    domain, verbose = True)
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
+                      'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 
+			'right': Bd})
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 240, finaltime = 10800): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topright')     
+
+for t in domain.evolve(yieldstep = 120, finaltime = 16200
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topright')     
+
+for t in domain.evolve(yieldstep = 60, finaltime = 21600
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topright')     
+    
+for t in domain.evolve(yieldstep = 120, finaltime = 27000
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topright')     
+
+for t in domain.evolve(yieldstep = 240, finaltime = 36000
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'topright')   
+  
+print 'That took %.2f seconds' %(time.time()-t0)
+
+print 'finished'
