Index: /anuga_work/production/perth/build_perth_250m.py
===================================================================
--- /anuga_work/production/perth/build_perth_250m.py	(revision 5782)
+++ /anuga_work/production/perth/build_perth_250m.py	(revision 5782)
@@ -0,0 +1,86 @@
+"""
+Script for building the elevation data to run tsunami inundation scenario
+for Perth, WA, Australia.
+
+Input: elevation data from project_250m.py
+Output: pts file stored in project_250m.topographies_dir 
+The run_perth.py is reliant on the output of this script.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.geospatial_data.geospatial_data import *
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+
+# Application specific imports
+import project_250m   # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+copy_code_files(project_250m.output_build_time_dir,__file__, 
+               dirname(project_250m.__file__)+sep+ project_250m.__name__+'.py' )
+
+start_screen_catcher(project_250m.output_build_time_dir)
+
+print 'USER:    ', project_250m.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
+#-------------------------------------------------------------------------------
+print"project_250m.poly_all",project_250m.poly_all
+print"project_250m.combined_dir_name",project_250m.combined_dir_name
+
+# input elevation directory filenames
+onshore_in_dir_name = project_250m.onshore_in_dir_name
+
+# output elevation directory filenames
+onshore_dir_name = project_250m.onshore_dir_name
+
+# 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)
+
+#creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
+
+#-------------------------------------------------------------------------------
+# Combine datasets into project_250m.combined_dir_name
+#-------------------------------------------------------------------------------
+
+print'create Geospatial data1 objects from topographies'
+G = Geospatial_data(file_name = onshore_dir_name + '.pts')
+
+print'clip combined geospatial object by bounding polygon'
+G_clipped = G.clip(project_250m.poly_all)
+
+print'export combined DEM file'
+if access(project_250m.topographies_dir,F_OK) == 0:
+    mkdir (project_250m.topographies_dir)
+G_clipped.export_points_file(project_250m.combined_dir_name + '.pts')
+#G_clipped.export_points_file(project_250m.combined_dir_name + '.txt') #Use for comparision in ARC
+
+import sys
+sys.exit()
Index: /anuga_work/production/perth/project_250m.py
===================================================================
--- /anuga_work/production/perth/project_250m.py	(revision 5782)
+++ /anuga_work/production/perth/project_250m.py	(revision 5782)
@@ -0,0 +1,249 @@
+"""Common filenames and locations for elevation, meshes and outputs.
+This script is the heart of all scripts in folder
+"""
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
+from anuga.utilities.system_tools import get_user_name, get_host_name
+from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+#------------------------------------------------------------------------------
+# Directory setup
+#------------------------------------------------------------------------------
+##Note: INUNDATIONHOME is the inundation directory, not the data directory.
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
+muxhome = getenv('MUXHOME')
+user = get_user_name()
+host = get_host_name()
+
+##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
+build_time = time+'_build'
+run_time = time+'_run'
+print 'gtime: ', gtime
+
+#------------------------------------------------------------------------------
+# Initial Conditions
+#------------------------------------------------------------------------------
+
+##Changed to reflect the modeled community (also sets up the directory system)
+state = 'western_australia'
+scenario_name = 'perth'
+scenario = 'perth_tsunami_scenario'
+##One or all can be changed each time the run_scenario script is excuted
+tide = 0 #0.6
+#event_number = 27255
+event_number = 27283
+alpha = 0.1
+friction=0.01
+starttime=0
+finaltime=80000 
+
+interior_mesh = 'all' # Can have 'all' or 'none' for Phase 2 study
+
+setup='final'  ## Can replace with trial or basic, this action will thin the mesh
+               ## so that the run script will take less time (hence less acurate) 
+
+if setup =='trial':
+    print'trial'
+    res_factor=10
+    time_thinning=48
+    yieldstep=240
+if setup =='basic': 
+    print'basic'
+    res_factor=4
+    time_thinning=12
+    yieldstep=120
+if setup =='final': 
+    print'final'
+    res_factor=1
+    time_thinning=4
+    yieldstep=60
+
+#------------------------------------------------------------------------------
+# Output Filename
+#------------------------------------------------------------------------------
+##Important to distiquish each run
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(event_number)+'_250m_' + str(interior_mesh) +'_'+str(user)
+
+#------------------------------------------------------------------------------
+# INPUT DATA
+#------------------------------------------------------------------------------
+
+##ELEVATION DATA## - used in build_perth.py
+## onshore data: format ascii grid with accompaning projection file
+grid_250m_2005 = 'grid' 
+
+##GAUGES## - used in get_timeseries.py
+gauge_name = 'perth.csv'
+gauge_name2 = 'thinned_MGA50.csv'
+
+##BOUNDING POLYGON## -used in build_boundary.py and run_perth.py respectively
+#NOTE: when files are put together must be in sequence - for ease go clockwise!
+#Check the run_perth.py for boundary_tags
+##thinned ordering file from Hazard Map: format index,latitude,longitude (with title)
+order_filename = 'thinned_boundary_ordering.txt'
+##landward bounding points
+landward = 'landward_bounding_polygon.txt'
+
+#------------------------------------------------------------------------------
+# OUTPUT ELEVATION DATA
+#------------------------------------------------------------------------------
+##Output filename for elevation
+## its a combination of all the data put together (utilisied in build_boundary)
+combined_name ='perth_combined_elevation_250m'
+
+#------------------------------------------------------------------------------
+# Directory Structure
+#------------------------------------------------------------------------------
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
+topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = anuga_dir+'topographies'+sep
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+boundaries_dir = anuga_dir+'boundaries'+ sep 
+output_dir = anuga_dir+'outputs'+sep
+gauges_dir = anuga_dir+'gauges'+sep
+meshes_dir = anuga_dir+'meshes'+sep
+
+#------------------------------------------------------------------------------
+# Location of input and output Data
+#------------------------------------------------------------------------------
+##where the input data sits
+onshore_in_dir_name = topographies_in_dir + grid_250m_2005
+
+##where the output data sits
+onshore_dir_name = topographies_dir + grid_250m_2005
+
+##where the combined file sits
+combined_dir_name = topographies_dir + combined_name
+
+##where the mesh sits (this is created during the run_perth.py)
+meshes_dir_name = meshes_dir + scenario_name + interior_mesh + '.msh'
+
+#where the boundary order 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_perth.py)
+landward_dir = boundaries_dir + landward
+
+#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
+output_build_time_dir = output_dir+build_time+dir_comment+sep #used for build_perth.py
+output_run_time_dir = output_dir+run_time+dir_comment+sep #used for run_perth.py
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#where the directory of the gauges sit
+gauges_dir_name = gauges_dir + gauge_name #used for get_timeseries.py
+gauges_dir_name2 = gauges_dir + gauge_name2 #used for get_timeseries.py
+
+#------------------------------------------------------------------------------
+# Interior region definitions
+#------------------------------------------------------------------------------
+
+#Land, to set the stage/water to be offcoast only
+poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
+
+#Initial bounding polygon for data clipping 
+poly_all = read_polygon(polygons_dir+'poly_all.csv')
+res_poly_all = 100000*res_factor
+
+#Area of Interest 1 (Fremantle)
+poly_aoi1 = read_polygon(polygons_dir+'CBD_coastal.csv')
+res_aoi1 = 500*res_factor
+
+#Area of Interest 2 (Rockingham)
+poly_aoi2 = read_polygon(polygons_dir+'rockingham_penguin.csv')
+res_aoi2 = 500*res_factor
+
+#Area of Interest 2 (garden Island)
+poly_aoi2a = read_polygon(polygons_dir+'garden.csv')
+res_aoi2a= 500*res_factor
+
+#Area of Interest 3 (geordie bay - record of tsunami impact)
+poly_aoi3 = read_polygon(polygons_dir+'geordie_bay.csv')
+res_aoi3 = 500*res_factor
+
+#Area of Interest 4 (sorrento - record of tsunami impact)
+poly_aoi4 = read_polygon(polygons_dir+'sorrento_gauge.csv')
+res_aoi4 = 500*res_factor
+
+#Area of Significance 1 (Garden Island and sand bank infront of Rockingham)
+poly_aos1 = read_polygon(polygons_dir+'garden_rockingham.csv')
+res_aos1 = 1000*res_factor
+
+#Area of Significance 2 (incorporate coastline of rottnest)
+poly_aos2 = read_polygon(polygons_dir+'rottnest_external.csv')
+res_aos2 = 1000*res_factor
+
+#Refined areas
+#Polygon designed to incorporate Dredge Area from Fremantle to
+#Rockingham the steep incline was making the mesh go to 0
+poly_aos3 = read_polygon(polygons_dir+'DredgeArea.csv')
+res_aos3 = 1000*res_factor
+
+#Shallow water 1
+poly_sw1 = read_polygon(polygons_dir+'internal_h20mORd3km.csv')
+res_sw1 = 25000*res_factor
+
+#Deep water (land of Rottnest, ANUGA does not do donuts!)
+poly_dw1 = read_polygon(polygons_dir+'rottnest_internal.csv')
+res_dw1 = 100000*res_factor
+
+# Combined all regions, must check that all are included!
+
+if interior_mesh =='all':
+    print'Mesh = all'
+    interior_regions = [[poly_aoi1,res_aoi1],[poly_aoi2,res_aoi2]
+                     ,[poly_aoi2a,res_aoi2a],[poly_aoi3,res_aoi3]
+                     ,[poly_aoi4,res_aoi4],[poly_aos1,res_aos1]
+                     ,[poly_aos2,res_aos2],[poly_aos3,res_aos3]
+                     ,[poly_sw1,res_sw1], [poly_dw1,res_dw1]]
+if interior_mesh =='none': 
+    print'Mesh = none'
+    interior_regions = []
+    
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+print 'min estimated number of triangles', trigs_min
+   
+
+#------------------------------------------------------------------------------
+# Clipping regions for export to asc and regions for clipping data
+#------------------------------------------------------------------------------
+
+#Geordie Bay extract ascii grid
+xminGeordie = 358000
+xmaxGeordie = 362000
+yminGeordie = 6458500
+ymaxGeordie = 6461000
+
+#Sorrento extract ascii grid
+xminSorrento = 379000
+xmaxSorrento = 382500
+yminSorrento = 6477000
+ymaxSorrento = 6480000
+
+#Fremantle extract ascii grid
+xminFremantle = 376000
+xmaxFremantle = 388000
+yminFremantle = 6449000
+ymaxFremantle = 6461000
+
+#Rockingham extract ascii grid
+xminRockingham = 373500
+xmaxRockingham = 385500
+yminRockingham = 6424000
+ymaxRockingham = 6433000
+
Index: /anuga_work/production/perth/run_perth_250m.py
===================================================================
--- /anuga_work/production/perth/run_perth_250m.py	(revision 5782)
+++ /anuga_work/production/perth/run_perth_250m.py	(revision 5782)
@@ -0,0 +1,238 @@
+"""Script for running tsunami inundation scenario for Perth, WA, Australia.
+
+The scenario is defined by a triangular mesh created from project_250m.polygon,
+the elevation data is combiled into a pts file through build_perth.py
+and a simulated tsunami is generated through an sts file from build_boundary.py
+
+Input: sts file (build_boundary.py for resepective event)
+       pts file (build_perth.py)
+       information from project_250m file
+Outputs: sww file stored in project_250m.output_run_time_dir 
+The export_results_all.py and get_timeseries.py is reliant
+on the outputs of this script
+
+Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
+Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# 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
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Field_boundary
+from Numeric import allclose
+from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
+from anuga_parallel.parallel_abstraction import get_processor_name
+from anuga.caching import myhash
+from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
+from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
+from polygon import Polygon_function
+    
+# Application specific imports
+import project_250m  # Definition of file names and polygons
+numprocs = 1
+myid = 0
+
+def run_model(**kwargs):
+    
+    #------------------------------------------------------------------------------
+    # Copy scripts to time stamped output directory and capture screen
+    # output to file
+    #------------------------------------------------------------------------------
+    print "Processor Name:",get_processor_name()
+
+    #copy script must be before screen_catcher
+
+    print 'output_dir',kwargs['output_dir']
+    
+    copy_code_files(kwargs['output_dir'],__file__, 
+             dirname(project_250m.__file__)+sep+ project_250m.__name__+'.py' )
+
+    store_parameters(**kwargs)
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+    
+    #-----------------------------------------------------------------------
+    # Domain definitions
+    #-----------------------------------------------------------------------
+
+    # Read in boundary from ordered sts file
+    urs_bounding_polygon=create_sts_boundary(os.path.join(project_250m.boundaries_dir_event,project_250m.scenario_name))
+
+    # Reading the landward defined points, this incorporates the original clipping
+    # polygon minus the 100m contour
+    landward_bounding_polygon = read_polygon(project_250m.landward_dir)
+
+    # Combine sts polyline with landward points
+    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
+    
+    # counting segments
+    N = len(urs_bounding_polygon)-1
+
+    #boundary tags refer to project_250m.landward 4 points equals 5 segments start at N
+    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
+
+    
+    #--------------------------------------------------------------------------
+    # Create the triangular mesh based on overall clipping polygon with a tagged
+    # boundary and interior regions defined in project_250m.py along with
+    # resolutions (maximal area of per triangle) for each polygon
+    #--------------------------------------------------------------------------
+
+    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
+    # causes problems with the ability to cache set quantity which takes alot of times
+        
+    print 'start create mesh from regions'
+
+    create_mesh_from_regions(bounding_polygon,
+                         boundary_tags=boundary_tags,
+                         maximum_triangle_area=project_250m.res_poly_all,
+                         interior_regions=project_250m.interior_regions,
+                         filename=project_250m.meshes_dir_name,
+                         use_cache=True,
+                         verbose=True)
+   
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    domain = Domain(project_250m.meshes_dir_name, use_cache=False, verbose=True)
+    print 'memory usage before del domain',mem_usage()
+       
+    print domain.statistics()
+    print 'triangles',len(domain)
+    
+    kwargs['act_num_trigs']=len(domain)
+
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+    print 'Setup initial conditions'
+
+    #following sets the stage/water to be offcoast only
+    IC = Polygon_function( [(project_250m.poly_mainland, 0)], default = kwargs['tide'],
+                             geo_reference = domain.geo_reference)
+    domain.set_quantity('stage', IC)
+    #domain.set_quantity('stage',kwargs['tide'] )
+    domain.set_quantity('friction', kwargs['friction']) 
+    
+    print 'Start Set quantity',kwargs['elevation_file']
+
+    domain.set_quantity('elevation', 
+                        filename = kwargs['elevation_file'],
+                        use_cache = False,
+                        verbose = True,
+                        alpha = kwargs['alpha'])
+    print 'Finished Set quantity'
+
+##   #------------------------------------------------------
+##    # Distribute domain to implement parallelism !!!
+##    #------------------------------------------------------
+##
+##    if numprocs > 1:
+##        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(kwargs['scenario_name'])
+    domain.set_datadir(kwargs['output_dir'])
+    domain.set_default_order(2) # Apply second order scheme
+    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
+    domain.set_store_vertices_uniquely(False)
+    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+    domain.tight_slope_limiters = 1
+    print 'domain id', id(domain)
+
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+   
+    boundary_urs_out=project_250m.boundaries_dir_event + sep + project_250m.scenario_name
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
+    
+    print 'Available boundary tags', domain.get_boundary_tags()
+    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
+                   domain, mean_stage= project_250m.tide,
+                   time_thinning=1,
+                   default_boundary=Bd,
+                   use_cache=True,
+                   verbose = True,
+                   boundary_polygon=bounding_polygon)
+
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': Bf}) 
+
+    kwargs['input_start_time']=domain.starttime
+
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = project_250m.yieldstep, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = False): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')
+
+    x, y = domain.get_maximum_inundation_location()
+    q = domain.get_maximum_inundation_elevation()
+
+    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+
+    print 'That took %.2f seconds' %(time.time()-t0)
+
+    #kwargs 'completed' must be added to write the final parameters to file
+    kwargs['completed']=str(time.time()-t0)
+     
+    store_parameters(**kwargs)
+     
+    print 'memory usage before del domain1',mem_usage()
+    
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    kwargs={}
+    kwargs['finaltime']=project_250m.finaltime
+    kwargs['output_dir']=project_250m.output_run_time_dir
+    kwargs['elevation_file']=project_250m.combined_dir_name+'.pts'
+    kwargs['scenario_name']=project_250m.scenario_name
+    kwargs['tide']=project_250m.tide
+    kwargs['alpha'] = project_250m.alpha
+    kwargs['friction']=project_250m.friction
+       
+
+    run_model(**kwargs)
+      
+   
