Index: anuga_work/production/broome/2007/export_results.py
===================================================================
--- anuga_work/production/broome/2007/export_results.py	(revision 5001)
+++ anuga_work/production/broome/2007/export_results.py	(revision 5001)
@@ -0,0 +1,79 @@
+import project_urs, os
+import sys
+
+from anuga.shallow_water.data_manager import sww2dem
+from os import sep
+
+time_dir = '20070615_063022_run_final_4.9_exmouth_nbartzis' # HAT exmouth
+#time_dir = '20070615_062730_run_final_0_exmouth_nbartzis' # MSL exmouth
+#time_dir = '20070618_064014_run_final_4.9_onslow_nbartzis' # HAT onslow
+#time_dir = '20070622_004008_run_final_0_onslow_nbartzis' # MSL onslow
+#time_dir = '20070615_063137_run_final_4.9_dampier_nbartzis' # HAT dampier 
+#time_dir = '20070621_002922_run_final_0_dampier_nbartzis' # MSL dampier 
+
+
+cellsize = 25
+timestep = None
+directory = project_urs.output_dir
+#name = directory + time_dir + sep + project.scenario_name
+name = directory+time_dir+sep+project_urs.scenario_name
+
+#is_parallel = True
+is_parallel = False
+if is_parallel == True: nodes = 4
+print 'output dir:', name
+
+var = [3]
+
+for which_var in var:
+    if which_var == 0:  # Stage
+        outname = name + '_stage'
+        quantityname = 'stage'
+
+    if which_var == 1:  # Absolute Momentum
+        outname = name + '_momentum_i1'
+        quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  
+
+    if which_var == 2:  # Depth
+        outname = name + '_depth'
+        quantityname = 'stage-elevation'  
+
+    if which_var == 3:  # Speed
+        outname = name + '_speed'
+        quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'  #Speed
+
+    if which_var == 4:  # Elevation
+        outname = name + '_elevation'
+        quantityname = 'elevation'  #Elevation
+
+    if is_parallel == True:
+    #    print 'is_parallel',is_parallel
+        for i in range(0,nodes):
+            namei = name + '_P%d_%d' %(i,nodes)
+            outnamei = outname + '_P%d_%d' %(i,nodes)
+            print 'start sww2dem for sww file %d' %(i)
+            sww2dem(namei, basename_out = outnamei,
+                        quantity = quantityname,
+                        timestep = timestep,
+                        cellsize = cellsize,      
+                        easting_min = project_urs.e_min_area,
+                        easting_max = project_urs.e_max_area,
+                        northing_min = project_urs.n_min_area,
+                        northing_max = project_urs.n_max_area,        
+                        reduction = max, 
+                        verbose = True,
+                        format = 'asc')
+    else:
+        print 'start sww2dem'
+        sww2dem(name, basename_out = outname,
+                    quantity = quantityname,
+                    timestep = timestep,
+                    cellsize = cellsize,      
+                    easting_min = project_urs.eastingmin,
+                    easting_max = project_urs.eastingmax,
+                    northing_min = project_urs.northingmin,
+                    northing_max = project_urs.northingmax,        
+                    reduction = max, 
+                    verbose = True,
+                    format = 'asc')
+
Index: anuga_work/production/broome/2007/project.py
===================================================================
--- anuga_work/production/broome/2007/project.py	(revision 5001)
+++ anuga_work/production/broome/2007/project.py	(revision 5001)
@@ -0,0 +1,181 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+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.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
+from anuga.utilities.system_tools import get_user_name
+
+# file and system info
+#---------------------------------
+codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+user = get_user_name()
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
+
+#time stuff
+time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
+#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
+
+#tide = -5.3
+tide = 0
+#tide = 4.9
+
+#Making assumptions about the location of scenario data
+state = 'western_australia'
+scenario_name = 'broome'
+scenario = 'broome_tsunami_scenario_2006'
+
+# onshore data provided by WA DLI
+onshore_name = 'Broome_topography_DLI' # original
+
+#island
+#island_name = 'rott_dli_ext' # original
+
+# offshore
+coast_name = 'Broome_coastline'
+
+offshore_name = 'Broome_Bathymetry'
+offshore_name1 = 'inferred_north'
+offshore_name2 = 'inferred_south'
+
+#final topo name
+combined_name ='broome_combined_elevation'
+combined_name1 ='broome_combined_elevation1'
+combined_name_unclipped1 ='broome_combined_elevation_unclipped1'
+combined_small_name = 'broome_combined_elevation_small'
+
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep+'Feb07'+sep
+topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
+topographies_time_dir = topographies_dir+build_time+sep
+
+#input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+#island_in_dir_name = topographies_in_dir + island_name
+
+coast_in_dir_name = topographies_in_dir + coast_name
+
+offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore_in_dir_name1 = topographies_in_dir + offshore_name1
+offshore_in_dir_name2 = topographies_in_dir + offshore_name2
+
+onshore_dir_name = topographies_dir + onshore_name
+
+#island_dir_name = topographies_dir + island_name
+
+coast_dir_name = topographies_dir + coast_name
+
+offshore_dir_name = topographies_dir + offshore_name
+offshore_dir_name1 = topographies_dir + offshore_name1
+offshore_dir_name2 = topographies_dir + offshore_name2
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_dir_name_unclipped1 = topographies_dir + combined_name_unclipped1
+combined_dir_name1 = topographies_dir + combined_name1
+combined_small_name_dir = topographies_dir + combined_small_name
+
+meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+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 = 'broome_3854_17042007'
+#boundaries locations
+boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep
+#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 + boundaries_source
+boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + boundaries_source
+#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
+
+#gauges
+beach_gauges = 'beach_gauges.csv'
+gauge_name = 'broome_gauges.csv'
+gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+beach_gauges = gauges_dir + 'beach_gauges.csv'
+gauges_dir_name = gauges_dir + gauge_name
+
+community_filename = gauges_dir + 'CHINS_v2.csv'
+community_broome = gauges_dir + 'community_broome.csv'
+
+buildings_filename = gauges_dir + 'Broome_res_Project.csv'
+buildings_filename_out = 'Broome_res_Project_modified.csv'
+
+###############################
+# Domain definitions
+###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+poly_all = read_polygon(polygons_dir+'extent_small.csv')
+#poly_all = read_polygon(polygons_dir+'extent.csv')
+res_factor = 2
+#res_poly_all = 500000
+res_poly_all = 150000*res_factor
+
+###############################
+# Interior region definitions
+###############################
+
+poly_0 = read_polygon(polygons_dir+'neg20_coast_contour_pts.csv')
+#res_0 = 100000
+res_0 = 20000*res_factor
+
+poly_1 = read_polygon(polygons_dir+'broome_north_coast_inside_extent.csv')
+#res_1 = 50000
+res_1 = 5000*res_factor
+
+poly_2 = read_polygon(polygons_dir+'broome_south_coast_inside_extent.csv')
+#res_2 = 50000
+res_2 = 5000*res_factor
+
+poly_3 = read_polygon(polygons_dir+'Broome_town_pts.csv')
+#res_3 = 20000
+res_3 = 2000*res_factor
+
+#poly_4 = read_polygon(polygons_dir+'Broome_inner_town_pts.csv')
+#res_4 = 5000
+#res_4 = 500*res_factor
+#assert zone == refzone
+
+#interior_regions = [[poly_0,res_0],[poly_1,res_1],[poly_2,res_2]
+                     #,[poly_3,res_3],[poly_4,res_4]]
+
+#trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+#print 'min number triangles', trigs_min
+
+#poly_mainland = read_polygon(polygons_dir+'Initial_Condition.csv')
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+e_min_area = 412900.0
+e_max_area = 422700.0
+n_min_area = 8008555.0
+n_max_area = 8021350.0
+
+
+
Index: anuga_work/production/broome/2007/run_broome.py
===================================================================
--- anuga_work/production/broome/2007/run_broome.py	(revision 5001)
+++ anuga_work/production/broome/2007/run_broome.py	(revision 5001)
@@ -0,0 +1,183 @@
+"""Script for running tsunami inundation scenario for Dampier, 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.output_run_time_dir 
+
+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 Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep, mkdir, access, F_OK
+from os.path import dirname, basename
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary, File_boundary, Reflective_boundary
+from anuga.shallow_water import Field_boundary
+from Numeric import allclose
+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
+from anuga.utilities.polygon import plot_polygons, polygon_area
+from anuga_parallel.parallel_abstraction import get_processor_name
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+if myid == 0:
+    copy_code_files(project.output_run_time_dir,__file__, 
+                 dirname(project.__file__)+sep+ project.__name__+'.py' )
+barrier()
+
+start_screen_catcher(project.output_run_time_dir, myid, numprocs)
+
+# filenames
+boundaries_name = project.scenario
+meshes_dir_name = project.meshes_dir_name+'.msh'
+boundaries_dir_name = project.boundaries_dir_name
+
+tide = project.tide
+
+# creates copy of code in output dir
+print "Processor Name:",get_processor_name()
+
+print 'USER: ', project.user
+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
+# tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#--------------------------------------------------------------------------
+
+if myid == 0:     
+
+    print 'start create mesh from regions'
+    create_mesh_from_regions(project.poly_all,
+#                         boundary_tags={'back': [1, 2], 'side': [0,3],
+#                                        'ocean': [4, 5, 6]},
+                         boundary_tags={'back': [1, 2], 'side': [0,3],
+                                        'ocean': [4]},
+                         maximum_triangle_area=project.res_poly_all,
+                         interior_regions=project.interior_regions,
+                         filename=meshes_dir_name,
+                         use_cache=False,
+                         verbose=True)
+
+# to sync all processors are ready 
+barrier()
+
+#-------------------------------------------------------------------------
+# Setup computational domain
+#-------------------------------------------------------------------------
+print 'Setup computational domain'
+#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
+domain = Domain(meshes_dir_name, use_cache=False, verbose=True)
+print domain.statistics()
+
+#-------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------
+if myid == 0:
+
+    print 'Setup initial conditions'
+
+    from polygon import Polygon_function
+    #following sets the stage/water to be offcoast only
+    IC = Polygon_function( [(project.poly_mainland, -1.0)], default = tide,
+                             geo_reference = domain.geo_reference)
+    domain.set_quantity('stage', IC)
+    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
+    print 'Start Set quantity'
+    print 'project.combined_dir_name_unclipped1',project.combined_dir_name_unclipped1+'.txt'
+    domain.set_quantity('elevation', 
+#                    filename = project.combined_dir_name+'.txt',
+                    filename = project.combined_dir_name_unclipped1+'.txt',
+#                    filename = project.combined_smaller_name_dir+'.xya',
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.1)
+    print 'Finished Set quantity'
+barrier()
+
+#------------------------------------------------------
+# Distribute domain to implement parallelism !!!
+#------------------------------------------------------
+
+if numprocs > 1:
+    domain=distribute(domain)
+
+#------------------------------------------------------
+# Set domain parameters
+#------------------------------------------------------
+
+domain.set_name(project.scenario_name)
+domain.set_datadir(project.output_run_time_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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+
+#-------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------
+print 'Available boundary tags', domain.get_boundary_tags()
+
+print 'Reading Boundary file'
+print 'domain id', id(domain)
+#Bf = File_boundary(boundaries_dir_name + '.sww',
+#                  domain, time_thinning=12, use_cache=True, verbose=True)
+Bf = Field_boundary(boundaries_dir_name + '.sww',
+                  domain, time_thinning=1, mean_stage=tide, 
+                  use_cache=True, verbose=True)
+
+print 'finished reading boundary file'
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+Bo = Dirichlet_boundary([tide+5.0,0,0])
+
+print'set_boundary'
+domain.set_boundary({'back': Br,
+                     'side': Bd,
+                     'ocean': Bf}) 
+#                     'ocean': Bd}) 
+print'finish set boundary'
+
+#----------------------------------------------------------------------------
+# Evolve system through time
+#----------------------------------------------------------------------------
+
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 60, finaltime = 30000):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'ocean')
+#    if allclose(t, 120):
+#        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
+
+#    if allclose(t, 1020):
+#        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
+
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
Index: anuga_work/production/broome/build_broome.py
===================================================================
--- anuga_work/production/broome/build_broome.py	(revision 5001)
+++ anuga_work/production/broome/build_broome.py	(revision 5001)
@@ -0,0 +1,163 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_urs.py
+The output sww file is stored in project_urs.output_time_dir 
+
+The scenario is defined by a triangular mesh created from project_urs.polygon,
+the elevation data and a simulated submarine landslide.
+
+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
+from anuga_parallel.parallel_abstraction import get_processor_name
+
+# Application specific imports
+import project_urs   # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+copy_code_files(project_urs.output_build_time_dir,__file__, 
+               dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
+
+start_screen_catcher(project_urs.output_build_time_dir)
+
+print 'USER:    ', project_urs.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_urs.combined_dir_name",project_urs.combined_dir_name
+
+# topography directory filenames
+onshore_in_dir_name = project_urs.onshore_in_dir_name
+onshore_in_dir_name1 = project_urs.onshore_in_dir_name1
+coast_in_dir_name = project_urs.coast_in_dir_name
+#island_in_dir_name = project_urs.island_in_dir_name
+offshore_in_dir_name = project_urs.offshore_in_dir_name
+offshore_in_dir_name1 = project_urs.offshore_in_dir_name1
+offshore_in_dir_name2 = project_urs.offshore_in_dir_name2
+offshore_in_dir_name3 = project_urs.offshore_in_dir_name3
+
+onshore_dir_name = project_urs.onshore_dir_name
+onshore_dir_name1 = project_urs.onshore_dir_name1
+coast_dir_name = project_urs.coast_dir_name
+#island_dir_name = project_urs.island_dir_name
+offshore_dir_name = project_urs.offshore_dir_name
+offshore_dir_name1 = project_urs.offshore_dir_name1
+offshore_dir_name2 = project_urs.offshore_dir_name2
+offshore_dir_name3 = project_urs.offshore_dir_name3
+
+# creates DEM from asc data
+print "creates DEMs from asc data"
+convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(onshore_in_dir_name1, basename_out=onshore_dir_name1, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(offshore_in_dir_name1, basename_out=offshore_dir_name1, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(offshore_in_dir_name2, basename_out=offshore_dir_name2, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(offshore_in_dir_name3, basename_out=offshore_dir_name3, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name,
+        use_cache=True, 
+        verbose=True)
+dem2pts(onshore_dir_name1, use_cache=True, verbose=True)
+#creates pts file for island DEM
+dem2pts(offshore_dir_name1, use_cache=True, verbose=True)
+dem2pts(offshore_dir_name2, use_cache=True, verbose=True)
+dem2pts(offshore_dir_name3, use_cache=True, verbose=True)
+
+print'create Geospatial data1 objects from topographies',onshore_dir_name + '.pts'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+print'create Geospatial data1a objects from topographies',onshore_dir_name1 + '.pts'
+G1a = Geospatial_data(file_name = onshore_dir_name1 + '.pts')
+print'create Geospatial data2 objects from coast', coast_in_dir_name + '.txt'
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt')
+#print'create Geospatial data3 objects from island'
+#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+print'create Geospatial data3 objects from offshore',offshore_in_dir_name + '.txt'
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt')
+print'create Geospatial data4 objects from offshore1',offshore_dir_name1 + '.pts'
+G_off1 = Geospatial_data(file_name = offshore_dir_name1 + '.pts')
+print'create Geospatial data4 objects from offshore2',offshore_dir_name2 + '.pts'
+G_off2 = Geospatial_data(file_name = offshore_dir_name2 + '.pts')
+print'create Geospatial data objects from offshore3',offshore_dir_name3 + '.pts'
+G_off3 = Geospatial_data(file_name = offshore_dir_name3 + '.pts')
+
+print'add all geospatial objects'
+G = G1 + G1a + G2 + G_off + G_off1 + G_off2 + G_off3
+
+print'clip combined geospatial object by bounding polygon'
+#G_clipped = G.clip(project_urs.poly_all)
+#FIXME: add a clip function to pts
+#print'shape of clipped data', G_clipped.get_data_points().shape
+
+print'export combined DEM file'
+if access(project_urs.topographies_dir,F_OK) == 0:
+    mkdir (project_urs.topographies_dir)
+print'export',project_urs.combined_dir_name+ '.txt'
+G.export_points_file(project_urs.combined_dir_name+ '.txt')
+
+'''
+print'split'
+G_small, G_other = G.split(.10,verbose=True)
+
+
+print 'export',project_urs.combined_dir_name + '.txt'
+G.export_points_file(project_urs.combined_dir_name + '.txt')
+print 'export', project_urs.combined_small_dir_name + '.txt'
+G_small.export_points_file(project_urs.combined_small_dir_name + '.txt')
+#G_clipped.export_points_file(project_urs.combined_dir_name + '.xya')
+
+
+
+print'export',project_urs.combined_dir_name+ '.txt'
+G_all=Geospatial_data(file_name = project_urs.combined_dir_name+ '.txt')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project_urs.combined_dir_name+'_small' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project_urs.combined_dir_name+'_other1' + '.xya')
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
+
+print 'boundaries_in_dir_name',project_urs.boundaries_in_dir_name
+
+urs_ungridded2sww(project_urs.boundaries_in_dir_name, project_urs.boundaries_in_dir_name,
+                  verbose=True, mint=4000, maxt=35000, zscale=1)
+
+
+'''
+
+
+
Index: anuga_work/production/broome/convert2latlon.py
===================================================================
--- anuga_work/production/broome/convert2latlon.py	(revision 5001)
+++ anuga_work/production/broome/convert2latlon.py	(revision 5001)
@@ -0,0 +1,43 @@
+"""
+Convert Easting Northing to lon/lat for ultimately get kml file
+"""
+
+def alter_file(filename,filenameout):
+
+    from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL
+    
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    
+    fid_out = open(filenameout, 'w')
+    s = '%s,%s,%s\n' %('lon','lat','stage')
+    fid_out.write(s)
+    for line in lines[1:]:
+       fields = line.split(',')
+       easting = float(fields[0])
+       northing = float(fields[1])
+       stage = float(fields[3])
+       lat, lon = UTMtoLL(northing,easting, 50)
+       s = '%f,%f,%f\n' %(lon,lat,stage)
+       fid_out.write(s)
+
+    fid_out.close() 
+
+    return
+
+import project
+from os import sep
+
+#timedir = '20070608_060316_run_final_1.5_onslow_nbartzis' # HAT 1 in Onslow 10000 yr
+timedir = '20070613_034754_run_final_1.5_dampier_nbartzis'
+#timedir = '20070906_061418_run_final_1.5_jsexton_exmouth_original'
+filename = project.outputdir + timedir + sep + 'stage0_interp.csv'
+#filename = project.outputdir + timedir + sep + 'stage5_interp.csv'
+#filename = project.outputdir+ timedir + sep + 'stage20_interp.csv'
+#filename = project.outputdir+  timedir + sep + 'stage80_interp.csv'
+filename_out = project.outputdir + timedir + sep + 'stage0_convert.csv'
+#filename_out = project.outputdir + timedir + sep + 'stage5_convert.csv'
+#filename_out = project.outputdir + timedir + sep + 'stage20_convert.csv'
+#filename_out = project.outputdir + timedir + sep + 'stage80_convert.csv'
+alter_file(filename,filename_out)
Index: anuga_work/production/broome/get_timeseries.py
===================================================================
--- anuga_work/production/broome/get_timeseries.py	(revision 5001)
+++ anuga_work/production/broome/get_timeseries.py	(revision 5001)
@@ -0,0 +1,81 @@
+"""
+Generate images of "gauges" for production run
+
+Inputs:
+
+production dirs: dictionary of production directories with a
+                 association to that simulation run, eg high tide,
+                 low tide and MSL.
+                    
+Outputs:
+
+* figures used for report stored in the report_figure directory
+
+"""
+from os import getcwd, sep, altsep, mkdir, access, F_OK
+import project
+from anuga.abstract_2d_finite_volumes.util import sww2timeseries
+
+# bulli event
+production_dirs = {'20070703_061553_run_final_0_dampier_nbartzis': ''}
+
+# 20070703_053716_run_final_4.9_exmouth_nbartzis  HAT exmouth
+# 20070703_053608_run_final_0_exmouth_nbartzis   # MSL exmouth
+# 20070703_053822_run_final_4.9_onslow_nbartzis  # HAT onslow
+# 20070703_061501_run_final_0_onslow_nbartzis  # MSL onslow
+# 20070703_062312_run_final_4.9_dampier_nbartzis  # HAT dampier 
+# 20070703_061553_run_final_0_dampier_nbartzis  # MSL dampier 
+
+is_parallel = True
+
+if is_parallel is True:
+    no_nodes = 4
+    
+if is_parallel is False:
+    
+    swwfiles = {}    
+    for label_id in production_dirs.keys():
+        
+        file_loc = project.output_dir + label_id + sep
+        swwfile = file_loc + project.scenario_name + '.sww'
+        swwfiles[swwfile] = label_id
+            
+        texname, elev_output = sww2timeseries(swwfiles,
+                                              project.beach_gauges,
+                                              production_dirs,
+                                              report = False,
+                                              plot_quantity = ['stage', 'speed'],
+                                              generate_fig = False,
+                                              surface = False,
+                                              time_min = None,
+                                              time_max = None,
+                                              time_unit = 'hours',
+                                              title_on = False,
+                                              use_cache = True,
+                                              verbose = True)
+
+if is_parallel is True:                                      
+
+    for count in range(no_nodes):
+        swwfiles = {} 
+        for label_id in production_dirs.keys():
+
+            file_loc = project.output_dir + label_id + sep
+            extra = '_P%i_%i' %(count,no_nodes)
+            swwfile = file_loc + project.scenario_name + extra + '.sww'
+            swwfiles[swwfile] = label_id
+
+            texname, elev_output = sww2timeseries(swwfiles,
+                                                  project.beach_gauges,
+                                                  production_dirs,
+                                                  report = False,
+                                                  plot_quantity = ['stage', 'speed'],
+                                                  generate_fig = False,
+                                                  surface = False,
+                                                  time_min = None,
+                                                  time_max = None,
+                                                  time_unit = 'hours',
+                                                  title_on = False,
+                                                  use_cache = True,
+                                                  verbose = True)
+                                      
Index: anuga_work/production/broome/locate_indigeneous.py
===================================================================
--- anuga_work/production/broome/locate_indigeneous.py	(revision 5001)
+++ anuga_work/production/broome/locate_indigeneous.py	(revision 5001)
@@ -0,0 +1,41 @@
+#
+# script to determine which communities are located in study region
+# 31 October 2006
+#
+
+def community_in_extent(filename,extent,zone,fid_out):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+
+    line1 = lines[0]
+    line11 = line1.split(',')
+    for i in range(len(line11)):
+        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
+        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
+        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
+
+    from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+    from anuga.utilities.polygon import plot_polygons, is_inside_polygon
+    
+    for line in lines[1:]:
+        fields = line.split(',')
+        point, zone1 = convert_from_latlon_to_utm(latitudes=[fields[lat_index]],\
+                                                 longitudes=[fields[lon_index]])
+        if zone1 == zone:
+            v = is_inside_polygon(point,extent,verbose=False)
+            if v == True:
+                p = point[0]
+                s = '%s, %.2f, %.2f\n' %(fields[name_index], p[0], p[1])
+                fid_out.write(s)
+    return
+
+import project
+
+filename = project.community_filename
+filename_out = project.community_broome
+fid_out = open(filename_out, 'w')
+extent = project.poly_all
+zone = 51
+
+community_in_extent(filename,extent,zone,fid_out)
Index: anuga_work/production/broome/plot_contour_data.py
===================================================================
--- anuga_work/production/broome/plot_contour_data.py	(revision 5001)
+++ anuga_work/production/broome/plot_contour_data.py	(revision 5001)
@@ -0,0 +1,65 @@
+"""Read in contour data file and plot max and min
+"""
+
+import project
+from pylab import plot, xlabel, ylabel, savefig, ion, close, axis, title, legend, grid
+from os import sep
+ 
+directory = project.outputdir
+
+timedir = '20070615_063022_run_final_4.9_exmouth_nbartzis' # HAT exmouth
+#timedir = '20070618_064014_run_final_4.9_onslow_nbartzis' # HAT onslow
+#timedir = '20070615_063137_run_final_4.9_dampier_nbartzis' # HAT dampier 
+
+def get_data(filename):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    stage = []
+    for line in lines[1:]:
+        fields = line.split(',')
+        stage.append(float(fields[2]))
+    return max(stage), min(stage)
+
+files = ['stage0_convert.csv','stage5_convert.csv','stage20_convert.csv','stage80_convert.csv']
+depth = [0,5,20,80]
+stage_max = []
+stage_min = []
+
+for thefile in files:
+
+    filename = directory + timedir + sep + thefile
+    max_stage, min_stage = get_data(filename)
+    stage_max.append(max_stage)
+    stage_min.append(min_stage)
+
+from anuga.abstract_2d_finite_volumes.util import greens_law
+from Numeric import arange
+d1 = 80.
+d2 = arange(d1,0.15,-0.1)
+h1 = max_stage
+#h1 = 4.0
+green = []
+for d in d2:
+    h2 = greens_law(d1,d,h1)
+    green.append(h2)
+    
+ion()
+plot(depth,stage_max,'ob',d2,green,'-g',depth,stage_max,'-b')
+xlabel('depth (m)') 
+ylabel('stage (m)')
+title('ANUGA modelled maximum stage versus Green\'s approximation \n \
+Magnitude 9.3 event (Java) at Broome')
+#title('ANUGA modelled maximum stage versus Green\'s approximation \n \
+#Magnitude 9.1 event (Sumba) at Broome')
+#title('ANUGA modelled maximum stage versus Green\'s approximation \n \
+#Magnitude 9.0 event (Sumba) at Broome')
+legend(['ANUGA','Green\'s law'])
+axis([-5,85,min(stage_max)*0.9,max(stage_max)*1.1])
+grid(True)
+#savefig('stratification_onslow_gun')
+savefig('stratification_exmouth_gun')
+#savefig('stratification_dampier_gun')
+
+close('all')
+
Index: anuga_work/production/broome/project_urs.py
===================================================================
--- anuga_work/production/broome/project_urs.py	(revision 5001)
+++ anuga_work/production/broome/project_urs.py	(revision 5001)
@@ -0,0 +1,216 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+Also includes origin for slump scenario.
+"""
+
+from os import sep, environ, getenv, getcwd,umask
+from os.path import expanduser, basename
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
+import sys
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
+from time import localtime, strftime, gmtime
+from anuga.utilities.system_tools import get_user_name, get_host_name
+
+codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir   
+user = get_user_name()
+host = get_host_name()
+#needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
+umask(002)
+
+
+#Making assumptions about the location of scenario data
+state = 'western_australia'
+scenario_name = 'broome'
+scenario = 'broome_tsunami_scenario_2006'
+
+#time stuff
+time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
+build_time = time+'_build'
+run_time = time+'_run'
+
+#tide = -5.3
+#tide = 0
+tide = 4.9
+
+#Maybe will try to make project a class to allow these parameters to be passed in.
+alpha = 0.1
+friction=0.01
+starttime=3600
+finaltime=25000
+setup='final'
+source='dampier'
+
+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
+
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
+
+# onshore data from 30m DTED level 2
+#onshore_name = 'Broome_topography_DLI' # original
+onshore_name = 'town_topo_10m' # original
+onshore_name1 = 'other_topo_250m' # original
+
+# offshore
+offshore_name = 'Broome_Bathymetry'
+#offshore_name1 = 'inferred_north'
+#offshore_name2 = 'inferred_south'
+offshore_name1 = 'north_250m'
+offshore_name2 = 'south_250m'
+offshore_name3 = 'cable_250m'
+coast_name = 'Broome_coastline'
+
+#final topo name
+combined_name ='broome_combined_elevation_z51_reduced'
+combined_name1 ='broome_combined_elevation1'
+combined_name_unclipped1 ='broome_combined_elevation_unclipped1'
+combined_small_name = 'broome_combined_elevation_small'
+
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
+
+topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep+'Feb07'+sep
+topographies_dir = anuga_dir+'topographies'+sep
+
+# input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+onshore_in_dir_name1 = topographies_in_dir + onshore_name1
+coast_in_dir_name = topographies_in_dir + coast_name
+offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore_in_dir_name1 = topographies_in_dir + offshore_name1
+offshore_in_dir_name2 = topographies_in_dir + offshore_name2
+offshore_in_dir_name3 = topographies_in_dir + offshore_name3
+
+onshore_dir_name = topographies_dir + onshore_name
+onshore_dir_name1 = topographies_dir + onshore_name1
+coast_dir_name = topographies_dir + coast_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
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_dir_name_unclipped1 = topographies_dir + combined_name_unclipped1
+combined_dir_name1 = topographies_dir + combined_name1
+combined_small_name_dir = topographies_dir + combined_small_name
+
+meshes_dir = anuga_dir+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+
+if source =='dampier':
+    boundaries_name = 'broome_3854_17042007' #Dampier gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep
+
+if source=='onslow':
+    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
+    
+if source=='exmouth':
+    boundaries_name = 'broome_3103_18052007' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
+    
+boundaries_in_dir_name = boundaries_in_dir + boundaries_name
+boundaries_dir = anuga_dir+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + boundaries_name
+
+#output locations
+output_dir = anuga_dir+'outputs'+sep
+output_build_time_dir = output_dir+build_time+'_'+source+sep
+#output_run_time_dir = output_dir +run_time+dir_comment+sep
+output_run_time_dir = output_dir +'20070703_062312_run_final_4.9_dampier_nbartzis'+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+gauge_name = 'broome_gauges.csv'
+gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+
+community_filename = gauges_dir + 'CHINS_v2.csv'
+community_broome = gauges_dir + 'community_broome.csv'
+
+buildings_filename = gauges_dir + 'Broome_res_Project.csv'
+buildings_filename_out = output_run_time_dir+'Broome_res_Project_modified.csv'
+
+###############################
+# Domain definitions
+###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+poly_all = read_polygon(polygons_dir+'extent_small.csv')
+# to test update of plot_polygons - now allows points and polygon to be
+# plotted together
+#point = read_polygon(polygons_dir+'neg20_coast_contour_pts.csv')
+#p2 = read_polygon(polygons_dir+'broome_north_coast_inside_extent.csv')
+#label=None
+#plot_polygons([poly_all,point,p2],['line','point','outside'],figname='test',label=label)
+#plot_polygons([poly_all,point],figname='test')
+#poly_all = read_polygon(polygons_dir+'extent.csv')
+res_poly_all = 150000*res_factor
+
+###############################
+# Interior region definitions
+###############################
+
+poly_0 = read_polygon(polygons_dir+'neg20_coast_contour_pts.csv')
+res_0 = 20000*res_factor
+
+poly_1 = read_polygon(polygons_dir+'broome_north_coast_inside_extent.csv')
+res_1 = 5000*res_factor
+
+#poly_2 = read_polygon(polygons_dir+'broome_south_coast_inside_extent.csv')
+poly_2 = read_polygon(polygons_dir+'broome_south_coast_not_town_pts_new.csv')
+res_2 = 5000*res_factor
+
+#poly_3 = read_polygon(polygons_dir+'Broome_town_pts.csv')
+poly_3 = read_polygon(polygons_dir+'broome_town_new_ptsa.csv')
+res_3 = 2000*res_factor
+
+#poly_4 = read_polygon(polygons_dir+'Broome_inner_town_pts.csv')
+poly_4 = read_polygon(polygons_dir+'broome_inner_town_new_ptsa.csv')
+res_4 = 500*res_factor
+
+interior_regions = [[poly_0,res_0],[poly_1,res_1],[poly_2,res_2]
+                     ,[poly_3,res_3],[poly_4,res_4]]
+                     
+boundary_tags={'back': [1, 2], 
+               'side': [0,3],
+               'ocean': [4]}
+
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
+
+poly_mainland = read_polygon(polygons_dir+'Initial_Condition.csv')
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 412000.0
+eastingmax = 423000.0
+northingmin = 8007000.0
+northingmax = 8022000.0
+
+# exporting asc grid for speed
+#eastingmin = 412900.0
+#eastingmax = 422700.0
+#northingmin = 8008555.0
+#northingmax = 8021350.0
Index: anuga_work/production/broome/run_boundary_points.py
===================================================================
--- anuga_work/production/broome/run_boundary_points.py	(revision 5001)
+++ anuga_work/production/broome/run_boundary_points.py	(revision 5001)
@@ -0,0 +1,35 @@
+"""Script for running tsunami inundation scenario for Dampier, 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.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+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
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+
+# # Related major packages
+# from anuga.shallow_water import Domain
+from anuga.shallow_water.data_manager import URS_points_needed_to_file, \
+     urs_ungridded2sww
+import project                 # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+poly = project.poly_all
+print "poly", poly
+print "project.boundaries_in_dir_name",project.boundaries_in_dir_name
+URS_points_needed_to_file(project.boundaries_in_dir_name,poly, export_csv=True, zone=51)
+#URS_points_needed_to_file(project.boundaries_dir_name5,poly, export_csv=True, zone=50)
+#URS_points_needed_to_file(project.boundaries_dir_name5,poly, export_csv=True)
+#urs_ungridded2sww(project.boundaries_in_dir_name2, project.boundaries_dir_name5, mean_stage= project.tide,  verbose=True)
Index: anuga_work/production/broome/run_broome_urs.py
===================================================================
--- anuga_work/production/broome/run_broome_urs.py	(revision 5001)
+++ anuga_work/production/broome/run_broome_urs.py	(revision 5001)
@@ -0,0 +1,298 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_urs.py
+The output sww file is stored in project_urs.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project_urs.polygon,
+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
+"""
+
+#------------------------------------------------------------------------------
+# 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 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
+
+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_api import distribute, numprocs, myid, barrier
+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 
+
+# Application specific imports
+import project_urs                 # Definition of file names and polygons
+
+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 kwargs
+    '''
+    print 'tide',kwargs['tide']
+    kwargs['est_num_trigs']=project_urs.trigs_min
+    kwargs['num_cpu']=numprocs
+    kwargs['host']=project_urs.host
+    kwargs['res_factor']=project_urs.res_factor
+    kwargs['starttime']=project_urs.starttime
+    kwargs['yieldstep']=project_urs.yieldstep
+    kwargs['finaltime']=project_urs.finaltime
+    
+    kwargs['output_dir']=project_urs.output_run_time_dir
+    kwargs['bathy_file']=project_urs.combined_dir_name+'.txt'
+#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
+    '''
+    print 'output_dir',kwargs['output_dir']
+    if myid == 0:
+        copy_code_files(kwargs['output_dir'],__file__, 
+                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
+
+        store_parameters(**kwargs)
+
+    barrier()
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+
+    # filenames
+#    meshes_dir_name = project_urs.meshes_dir_name+'.msh'
+
+    # creates copy of code in output dir
+    print 'min triangles', project_urs.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
+    # tagged
+    # boundary and interior regions defined in project_urs.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
+    if myid == 0:
+    
+        print 'start create mesh from regions'
+
+        create_mesh_from_regions(project_urs.poly_all,
+                         boundary_tags=project_urs.boundary_tags,
+                             maximum_triangle_area=project_urs.res_poly_all,
+                             interior_regions=project_urs.interior_regions,
+                             filename=project_urs.meshes_dir_name+'.msh',
+                             use_cache=False,
+                             verbose=True)
+    barrier()
+    
+
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
+    #above don't work
+    domain = Domain(project_urs.meshes_dir_name+'.msh', 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
+    #-------------------------------------------------------------------------
+    if myid == 0:
+
+        print 'Setup initial conditions'
+
+        from polygon import Polygon_function
+        #following sets the stage/water to be offcoast only
+#        IC = Polygon_function( [(project_urs.poly_mainland, -1.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['bathy_file']
+
+        domain.set_quantity('elevation', 
+                            filename = kwargs['bathy_file'],
+                            use_cache = True,
+                            verbose = True,
+                            alpha = kwargs['alpha'])
+        print 'Finished Set quantity'
+    barrier()
+
+
+    #------------------------------------------------------
+    # Distribute domain to implement parallelism !!!
+    #------------------------------------------------------
+
+    if numprocs > 1:
+        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(kwargs['aa_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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+    #domain.tight_slope_limiters = 1
+
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
+
+    Bf = Field_boundary(kwargs['boundary_file'],
+                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
+                    use_cache=True, verbose=True)
+                    
+    kwargs['input_start_time']=domain.starttime
+
+    print 'finished reading boundary file'
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
+
+    print'set_boundary'
+
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': Bf}) 
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+
+    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
+                       ,skip_initial_step = True): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+    
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = True): 
+        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)
+    
+    if myid==0:
+        store_parameters(**kwargs)
+    barrier
+    
+    print 'memory usage before del domain1',mem_usage()
+    
+def export_model(**kwargs):    
+    #store_parameters(**kwargs)
+    
+    print 'memory usage before del domain',mem_usage()
+    #del domain
+    print 'memory usage after del domain',mem_usage()
+    
+    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
+    print'swwfile',swwfile
+    
+    export_grid(swwfile, extra_name_out = 'town',
+            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            timestep = None,
+            reduction = max,
+            cellsize = 25,
+            NODATA_value = -9999,
+            easting_min = project_urs.eastingmin,
+            easting_max = project_urs.eastingmax,
+            northing_min = project_urs.northingmin,
+            northing_max = project_urs.northingmax,
+            verbose = False,
+            origin = None,
+            datum = 'GDA94',
+            format = 'asc')
+            
+    inundation_damage(swwfile+'.sww', project_urs.buildings_filename, 
+                      project_urs.buildings_filename_out)
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    kwargs={}
+    kwargs['est_num_trigs']=project_urs.trigs_min
+    kwargs['num_cpu']=numprocs
+    kwargs['host']=project_urs.host
+    kwargs['res_factor']=project_urs.res_factor
+    kwargs['starttime']=project_urs.starttime
+    kwargs['yieldstep']=project_urs.yieldstep
+    kwargs['finaltime']=project_urs.finaltime
+    
+    kwargs['output_dir']=project_urs.output_run_time_dir
+    kwargs['bathy_file']=project_urs.combined_dir_name+'.txt'
+#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
+    kwargs['file_name']=project_urs.home+'detail.csv'
+    kwargs['aa_scenario_name']=project_urs.scenario_name
+    kwargs['ab_time']=project_urs.time
+    kwargs['res_factor']= project_urs.res_factor
+    kwargs['tide']=project_urs.tide
+    kwargs['user']=project_urs.user
+    kwargs['alpha'] = project_urs.alpha
+    kwargs['friction']=project_urs.friction
+    kwargs['time_thinning'] = project_urs.time_thinning
+    kwargs['dir_comment']=project_urs.dir_comment
+
+    run_model(**kwargs)
+      
+    if myid==0:
+        export_model(**kwargs)
+    barrier
Index: anuga_work/production/broome/run_building_inundation.py
===================================================================
--- anuga_work/production/broome/run_building_inundation.py	(revision 5001)
+++ anuga_work/production/broome/run_building_inundation.py	(revision 5001)
@@ -0,0 +1,42 @@
+""" Script for calculating maximum inundation at building locations
+
+    Inputs:
+
+    swwfile             - name of sww file
+                        - assume that all conserved quantities have been stored
+                    
+    buildings_filename  - name of file containing building data, sourced
+                          from NBED
+
+    Output:
+    
+    - building_filename (.csv) augmented to include maximum depth, momentum and
+      velocity and stored in same location as swwfile
+      Name = augmented_buildings.csv
+      
+    """
+
+from os import sep
+
+from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
+import project_urs
+
+# Inputs
+#time_dir = '20070615_063022_run_final_4.9_exmouth_nbartzis' # HAT exmouth
+#time_dir = '20070615_062730_run_final_0_exmouth_nbartzis' # MSL exmouth
+#time_dir = '20070618_064014_run_final_4.9_onslow_nbartzis' # HAT onslow
+#time_dir = '20070622_004008_run_final_0_onslow_nbartzis' # MSL onslow
+time_dir = '20070615_063137_run_final_4.9_dampier_nbartzis' # HAT dampier 
+#time_dir = '20070621_002922_run_final_0_dampier_nbartzis' # MSL dampier
+
+file_loc = project_urs.output_dir + time_dir + sep 
+swwfile = file_loc + project_urs.scenario_name #+ '.sww'
+buildings_filename = project_urs.buildings_filename
+buildings_filename_out = file_loc + project_urs.buildings_filename_out
+
+#add_depth_and_momentum2csv(swwfile, buildings_filename, buildings_filename_out)
+inundation_damage(swwfile, buildings_filename, buildings_filename_out)
+
+
+print '\n Augmented building file written to %s \n' \
+      %buildings_filename_out
