Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py	(revision 7628)
@@ -0,0 +1,96 @@
+"""Build the elevation data to run a storm surge inundation scenario 
+for Mandurah, Western Australia, Australia.
+
+Input: elevation data from project.py
+Output: pts file stored in project.topographies_dir 
+The run_model.py is reliant on the output of this script.
+
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+from os.path import join
+
+# Related major packages
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
+from anuga.shallow_water.data_manager import dem2pts
+from anuga.shallow_water.data_manager import start_screen_catcher
+from anuga.shallow_water.data_manager import copy_code_files
+from anuga.geospatial_data.geospatial_data import Geospatial_data
+
+# Application specific imports
+from setup_model_GEMS import project_GEMS   # Definition of file names and polygons
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+copy_code_files(project_GEMS.output_build,__file__, 
+				os.path.dirname(project_GEMS.__file__)+os.sep+\
+				project_GEMS.__name__+'.py' )
+start_screen_catcher(project_GEMS.output_build)
+
+#------------------------------------------------------------------------------
+# 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.bounding_polygon', project_GEMS.bounding_polygon
+print 'project.combined_elevation_basename', project_GEMS.combined_elevation_basename
+
+# Create Geospatial data from ASCII files
+geospatial_data = {}
+if not project_GEMS.ascii_grid_filenames == []:
+    for filename in project_GEMS.ascii_grid_filenames:
+        absolute_filename = join(project_GEMS.topographies_folder, filename)
+        convert_dem_from_ascii2netcdf(absolute_filename,
+                                      basename_out=absolute_filename,
+                                      use_cache=True,  
+                                      verbose=True)
+        dem2pts(absolute_filename, use_cache=True, verbose=True)
+
+        G_grid = Geospatial_data(file_name=absolute_filename+'.pts',
+                                                    verbose=True)
+        print 'Clip geospatial object'
+        geospatial_data[filename] = G_grid.clip(project_GEMS.bounding_polygon)
+
+# Create Geospatial data from TXT files
+if not project_GEMS.point_filenames == []:
+    for filename in project_GEMS.point_filenames:
+        absolute_filename = join(project_GEMS.topographies_folder, filename)
+        G_points = Geospatial_data(file_name=absolute_filename,
+                                                    verbose=True)
+        print 'Clip geospatial object'
+        geospatial_data[filename] = G_points.clip(project_GEMS.bounding_polygon)
+
+print 'Geospatial data objects: ', geospatial_data.keys()
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+
+print 'Add geospatial objects' # except', project_GEMS.offshore_name5
+
+print 'project_GEMS.elevation_clip_box', project_GEMS.elevation_clip_box
+
+G1 = geospatial_data['Man_25m'].clip(project_GEMS.elevation_clip_box)
+G2 = geospatial_data['man10m_ss'].clip_outside(project_GEMS.elevation_clip_box)
+
+G = G1 + G2
+	
+print 'Export combined DEM file'
+G.export_points_file(project_GEMS.combined_elevation + '.pts')
+
+print 'Do txt version too'
+# Use for comparision in ARC
+G.export_points_file(project_GEMS.combined_elevation + '.txt')
+G1.export_points_file(join(project_GEMS.topographies_folder, (str(project_GEMS.ascii_grid_filenames[0]) +'_export.txt')))
+G2.export_points_file(join(project_GEMS.topographies_folder, (str(project_GEMS.ascii_grid_filenames[1]) + '_export.txt')))
Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/file_length.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/file_length.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/file_length.py	(revision 7628)
@@ -0,0 +1,22 @@
+"""Function to return the length of a file.
+
+Intended to be used for simplfying the creation of boundary_tags in
+run_model.py
+
+Input: filename
+Returns: number of lines in file
+"""
+
+def file_length(in_file):
+    '''Function to return the number of lines in a file.
+
+    in_file: Path to the file to get number of lines in.
+    
+    Returns: number of lines in file
+    '''
+    
+    fid = open(in_file)
+    data = fid.readlines()
+    fid.close()
+    return len(data)
+
Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/get_runup.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/get_runup.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/get_runup.py	(revision 7628)
@@ -0,0 +1,69 @@
+# Standard modules
+from os import sep
+import os
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from os.path import join
+from anuga.utilities.polygon import read_polygon
+from shutil import copy
+import time
+import sys
+# Related major packages
+from anuga.shallow_water.data_manager import get_maximum_inundation_data
+from anuga.shallow_water.data_manager import csv2polygons
+
+# Application specific imports
+from setup_model import project  # Definition of file names and polygons
+import csv
+directory = project.output_folder
+
+time_dir1 = '20090508_150215_run_final_0_51469_lfountai'
+time_dir2 = '20090511_161526_run_final_1.1_51469_kvanputt'
+time_dir3 = '20090511_165539_run_final_0_50863_lfountai'
+time_dir4 = '20090521_220101_run_final_1.1_50863_kvanputt'
+time_dir5 = '20090518_154710_run_final_0_50994_lfountai'
+time_dir6 = '20090519_160510_run_final_1.1_50994_lfountai'
+time_dir7 = '20090522_164526_run_final_0_51392_lfountai'
+time_dir8 = '20090522_164640_run_final_1.1_51392_lfountai'
+time_dir9 = '20090522_164948_run_final_0_51423_lfountai'
+time_dir10 = '20090522_165600_run_final_1.1_51423_lfountai'
+
+time_dirs = [time_dir1, time_dir2, time_dir3, time_dir4, time_dir5,
+             time_dir6, time_dir7, time_dir8, time_dir9, time_dir10]
+
+
+# polygons is a dictionary of polygon data indexed by 'id'
+# values is a dictionary of polygon names indexed by 'id'
+polygons, values = csv2polygons(project.images)
+print 'Creating %d highest runup' % len(polygons)
+
+# create a dictionary {<polygon_name>: <polygon_data>, ...}
+area_dictionary = {}
+for id in polygons:
+    area_dictionary[values[id]] = polygons[id] 
+
+# for each time directory ...
+for time_dir in time_dirs:
+    print 'time_dir ',time_dir
+    basename = join(directory, time_dir, project.scenario_name)
+
+    out_file = basename+'_runup.csv'
+    output_file = open(out_file, 'w')
+    s = 'area,q,x,y\n' 
+    output_file.write(s)
+
+    # iterate over polygon names
+    for which_area in values.itervalues():
+        print 'which_area = %s' % which_area
+        polygon = area_dictionary[which_area]
+        
+        q, loc = get_maximum_inundation_data(basename+'.sww',
+                                             polygon=polygon, verbose=True)
+
+        print '%s, %.6f, %.6f, %.6f\n' %(str(which_area), q, loc[0],loc[1])
+        s = '%s, %.6f, %.6f, %.6f\n' %(str(which_area), q, loc[0],loc[1])
+
+        output_file.write(s)
+
+        print 'output file printed to: %s' %out_file
+    output_file.close()
Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/get_timeseries.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/get_timeseries.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/get_timeseries.py	(revision 7628)
@@ -0,0 +1,31 @@
+"""
+Generate time series of nominated "gauges".
+Input: sww file from run_perth.py
+       gauges project.gauge_filename
+Output: csv files stage,speed,depth,elevation over time
+Stored in the 'outputs_dir' folder for respective .sww file
+
+Note:
+Can run different sww files at the same time
+make sure if there is a second sww file that it is placed into a folder called sww2
+Can run different gauges at the same time - ie testing boundary index point
+"""
+
+from os import sep, rename, listdir, system
+from os.path import join
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
+import project
+
+directory = project.output_folder
+
+time_dir1 = '20091229_100831_run_storm_surge_final_0_alby_waves_lfountai'
+
+time_dirs = [time_dir1]
+
+for time_dir in time_dirs:
+    name = join(directory,time_dir,project.scenario_name)
+    gauge = project.gauges
+    sww2csv_gauges(name+'.sww',gauge,
+                   quantities = ['stage','speed','depth','elevation'],
+                   verbose=True)
+
Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py	(revision 7628)
@@ -0,0 +1,259 @@
+"""
+This file contains all your file and directory definitions 
+for elevation, meshes and outputs.
+"""
+
+import os
+from anuga.utilities.system_tools import get_user_name, get_host_name
+from time import localtime, strftime, gmtime
+from os.path import join, exists
+
+
+#-------------------------------------------------------------------------------
+# Directory setup
+#-------------------------------------------------------------------------------
+
+# this section needs to be updated to reflect the modelled community.
+# Note, the user needs to set up the directory system accordingly
+state = 'western_australia'
+#scenario_name = 'case_a'
+scenario_name = 'alby_waves'
+scenario_folder = 'mandurah_storm_surge_scenario_2009'
+
+#-------------------------------------------------------------------------------
+# Initial Conditions
+#-------------------------------------------------------------------------------
+
+# Model specific parameters.
+# One or all can be changed each time the run_model script is executed
+
+
+central_meridian = None # Central meridian for projection (optional)
+zone = 50               # UTM zone for projection
+
+event = 'alby_waves'
+#event = 'case_a'        # case_a - baseline historic TC Alby event
+                        # case_b - worst case scenarion current climate (current coastline)
+                        # case_c - worst case 1.1m SLR current climate (current coastline)
+                        # case_d - worst case 0.8m SLR current climate (current coastline)
+                        # case_e - worst case 0.5m SLR current climate (current coastline)
+tide = 0                # Mean Sea Level = 0,
+						# Highest Astronomical Tide = 0.6 m for Bunbury - CHECK THIS! 
+alpha = 0.1             # smoothing parameter for mesh
+friction = 0.01           # manning's friction coefficient
+starttime = 0             # start time for simulation
+finaltime = 86400         # final time for simulation - 24 hours for TC Alby
+
+setup = 'storm_surge_final'         # This can be one of four values
+						#    trial - coarsest mesh, fast
+						#    basic - coarse mesh
+						#    final - fine mesh, slowest
+						#    storm_surge_final - as 'final' but with a longer yieldstep (12 mins)
+
+#-------------------------------------------------------------------------------
+# Output filename
+#
+# Your output filename should be unique between different runs on different data.
+# The list of items below will be used to create a file in your output directory.
+# Your user name and time+date will be automatically added.  For example,
+#     [setup, tide, event_number]
+# will result in a filename like
+#     20090212_091046_run_final_0_27283_rwilson
+#-------------------------------------------------------------------------------
+
+output_comment = [setup, tide, event, '_GEMS_compare']    # event_number will have to
+												# change to something relevent
+												# for storm surge
+
+#-------------------------------------------------------------------------------
+# Input Data
+#-------------------------------------------------------------------------------
+
+# ELEVATION DATA
+# Used in build_elevation.py
+# Format for ascii grids, as produced in ArcGIS + a projection file
+##ascii_grid_filenames = ['m_peel.asc',
+##                        'm_harvey.asc',
+##                        'm_peel_10m.asc'
+##                        'm_harvey_10m.asc']
+ascii_grid_filenames = ['Man_25m',		# this is the latest 25m DEM from GEMS
+						'man10m_ss']	# this is to fill in areas not covered by the 
+                                        # GEMS grid and is derived from the 10m resampled LiDAR from DoT,
+                                        # clipped to the bounding polygon
+
+# Format for point is x,y,elevation (with header)
+point_filenames = []
+##			['MA-46893-SNDS_AHD.csv',
+##                   'MS0205HY_AHD.csv',
+##                   'MS0404_AHD.csv',
+##                   'YU0403HY_AHD.csv']
+
+# BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
+# Used in build_elevation.py & run_model.py
+# Format for points easting,northing (no header)
+bounding_polygon_filename = 'bounding_polygon_ss.csv'
+bounding_polygon_maxarea = 50000
+
+# INTERIOR REGIONS -  for designing the mesh
+# Used in run_model.py
+# Format for points easting,northing (no header)                    
+interior_regions_data = [['intermediate.csv', 2500],
+						['area_of_interest.csv', 100],
+                        ['canals_nth_2.csv', 100],  # Canal regions used to force mesh to fit canals
+                        ['canals_nth_3.csv', 100],   
+                        ['canals_nth_4.csv', 100], 
+                        ['canals_nth_5.csv', 100],
+                        ['canals_nth_6.csv', 100], 
+                        ['canals_nth_7.csv', 100],
+                        ['canals_nth_8.csv', 100], 
+                        ['canals_nth_9.csv', 100],
+                        ['canals_nth_10.csv', 100], 
+                        ['canals_nth_11.csv', 100],
+                        ['canals_nth_12.csv', 100], 
+                        ['canals_sth_1.csv', 100],
+                        ['canals_sth_2.csv', 100], 
+                        ['canals_sth_3.csv', 100],
+                        ['canals_sth_4.csv', 100], 
+                        ['canals_sth_5.csv', 100],
+                        ['canals_sth_6.csv', 100]]
+                        
+# LAND - used to set the initial stage/water to be offcoast only
+# Used in run_model.py.  Format for points easting,northing (no header)
+land_initial_conditions_filename = [['initial_conditions.csv', 0]]
+
+# GEMS order filename
+# Format is index,northing, easting, elevation (without header)
+gems_order_filename = 'gems_boundary_order_thinned.csv'
+
+# GAUGES - for creating timeseries at a specific point
+# Used in get_timeseries.py.  
+# Format easting,northing,name,elevation (with header)
+gauges_filename = 'gauges_mandurah_edited.csv'
+
+# BUILDINGS EXPOSURE - for identifying inundated houses
+# Used in run_building_inundation.py
+# Format latitude,longitude etc (geographic)
+building_exposure_filename = '' # from NEXIS
+
+# Landward bounding points
+# Format easting,northing (no header)
+landward_boundary_filename = 'landward_boundary.csv'
+
+### MUX input filename.
+### If a meta-file from EventSelection is used, set 'multi-mux' to True.
+### If a single MUX stem filename (*.grd) is used, set 'multi-mux' to False.
+####mux_input_filename = event_number # to be found in event_folder
+##                                    # (ie boundaries/event_number/)
+####multi_mux = False
+##mux_input_filename = 'event.list'
+##multi_mux = True
+
+#-------------------------------------------------------------------------------
+# Clipping regions for export to asc and regions for clipping data
+# Final inundation maps should only be created in regions of the finest mesh
+#-------------------------------------------------------------------------------
+
+# Elevation clip box - used when there are multiple overlapping datasets and you want to use 
+# one over another in an area of overlap
+elevation_clip_box_filename = 'elevation_clip_box.csv'
+
+################################################################################
+################################################################################
+####         NOTE: NOTHING WOULD NORMALLY CHANGE BELOW THIS POINT.          ####
+################################################################################
+################################################################################
+
+# Get system user and host names.
+# These values can be used to distinguish between two similar runs by two
+# different users or runs by the same user on two different machines.
+user = get_user_name()
+host = get_host_name()
+
+# Environment variable names.
+# The inundation directory, not the data directory.
+ENV_INUNDATIONHOME = 'INUNDATIONHOME'
+
+#-------------------------------------------------------------------------------
+# Output Elevation Data
+#-------------------------------------------------------------------------------
+
+# Output filename for elevation
+# this is a combination of all the data generated in build_elevation.py
+combined_elevation_basename = scenario_name + 'GEMS_combined_elevation'
+
+#-------------------------------------------------------------------------------
+# Directory Structure
+#-------------------------------------------------------------------------------
+
+# determines time for setting up output directories
+time = strftime('%Y%m%d_%H%M%S', localtime()) 
+gtime = strftime('%Y%m%d_%H%M%S', gmtime()) 
+build_time = time + '_build'
+run_time = time + '_run_'
+
+# create paths generated from environment variables.
+home = join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
+	
+# check various directories/files that must exist
+anuga_folder = join(home, state, scenario_folder, 'anuga')
+topographies_folder = join(anuga_folder, 'topographies')
+polygons_folder = join(anuga_folder, 'polygons')
+boundaries_folder = join(anuga_folder, 'boundaries')
+output_folder = join(anuga_folder, 'outputs')
+gauges_folder = join(anuga_folder, 'gauges')
+event_folder = join(boundaries_folder, str(event))
+
+#-------------------------------------------------------------------------------
+# Location of input and output data
+#-------------------------------------------------------------------------------
+
+# Convert the user output_comment to a string for run_model.py
+output_comment = ('_'.join([str(x) for x in output_comment if x != user])
+				  + '_' + user)
+
+# The absolute pathname of the all elevation, generated in build_elevation.py
+combined_elevation = join(topographies_folder, combined_elevation_basename)
+
+
+# The pathname for the GEMS order points
+if gems_order_filename:
+	gems_order = join(boundaries_folder, gems_order_filename)
+
+# The absolute pathname for the landward points of the bounding polygon,
+# Used within run_model.py)
+if landward_boundary_filename:
+	landward_boundary = join(boundaries_folder, landward_boundary_filename)
+
+# The absolute pathname for the .sts file, generated in build_boundary.py
+event_sts = join(event_folder, scenario_name)
+
+# The absolute pathname for the output folder names
+# Used for build_elevation.py
+output_build = join(output_folder, build_time) + '_' + str(user) 
+# Used for run_model.py
+output_run = join(output_folder, run_time) + output_comment 
+# Used by post processing
+output_run_time = join(output_run, scenario_name) 
+
+# The absolute pathname of the mesh, generated in run_model.py
+meshes = join(output_run, scenario_name) + 'GEMS.msh'
+
+# The absolute pathname for the gauges file
+# Used for get_timeseries.py
+if gauges_filename:
+	gauges = join(gauges_folder, gauges_filename)       
+
+# The absolute pathname for the building file
+# Used for run_building_inundation.py
+if building_exposure_filename:
+	building_exposure = join(gauges_folder, building_exposure_filename)
+
+### The absolute pathname for the image file
+### Used for get_runup.py
+##if images_filename:
+##    images = join(polygons_folder, images_filename)
+
+### full path to where MUX files (or meta-files) live
+##mux_input = join(event_folder, mux_input_filename)
+
+
Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py	(revision 7628)
@@ -0,0 +1,171 @@
+"""Run a storm surge inundation scenario for Mandurah, WA, Australia using
+input from the third party consultant GEMS.
+
+The scenario is defined by a triangular mesh created from project.polygon, the
+elevation data is compiled into a pts file through build_elevation.py and a
+simulated storm surge is generated from GEMS data.
+
+Input: sts files (csv format time series at ocean boundary for respective event)
+       pts file (build_elevation.py)
+       information from project file
+Outputs: sww file stored in project.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
+Nariman Habili, Leharne Fountain - 2009
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+##import os.path
+import time
+##from time import localtime, strftime, gmtime
+
+# Related major packages
+##from Scientific.IO.NetCDF import NetCDFFile
+import numpy as num
+
+from anuga.interface import create_domain_from_regions
+from anuga.interface import Transmissive_stage_zero_momentum_boundary
+from anuga.interface import Dirichlet_boundary
+from anuga.interface import Reflective_boundary
+from anuga.interface import Field_boundary
+from anuga.interface import create_sts_boundary
+##from anuga.interface import csv2building_polygons
+from file_length import file_length
+
+from anuga.shallow_water.data_manager import start_screen_catcher
+from anuga.shallow_water.data_manager import copy_code_files
+##from anuga.shallow_water.data_manager import urs2sts
+from anuga.utilities.polygon import read_polygon, Polygon_function
+
+# Application specific imports
+from setup_model_GEMS import project_GEMS
+##import build_urs_boundary as bub
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file. Copy script must be before screen_catcher
+#-------------------------------------------------------------------------------
+
+copy_code_files(project_GEMS.output_run, __file__, 
+                os.path.join(os.path.dirname(project_GEMS.__file__),
+                             project_GEMS.__name__+'.py'))
+start_screen_catcher(project_GEMS.output_run, 0, 1)
+
+#-------------------------------------------------------------------------------
+# Create the computational domain based on overall clipping polygon with
+# a tagged boundary and interior regions defined in project_GEMS.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+print 'Create computational domain'
+
+### Create the STS file
+##print 'project.mux_data_folder=%s' % project.mux_data_folder
+##if not os.path.exists(project.event_sts + '.sts'):
+##    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
+
+### Read in boundary from ordered sts file
+##gems_boundary = create_sts_boundary(project_GEMS.event_sts)
+
+#reading the GEMS boundary points
+gems_boundary = read_polygon(project_GEMS.gems_order)
+##gems_boundary = create_sts_boundary(project_GEMS.event_sts)
+
+# Reading the landward defined points of the original bounding polygon
+landward_boundary = read_polygon(project_GEMS.landward_boundary)
+
+# Combine GEMS input boundary polyline with landward points
+bounding_polygon_gems = gems_boundary + landward_boundary
+
+# Number of boundary segments
+num_ocean_segments = len(gems_boundary) - 1
+# Number of landward_boundary points
+num_land_points = file_length(project_GEMS.landward_boundary)
+
+# Boundary tags refer to project_GEMS.landward_boundary
+# 4 points equals 5 segments start at N
+boundary_tags={'back': range(num_ocean_segments+1,
+                             num_ocean_segments+num_land_points),
+               'side': [num_ocean_segments,
+                        num_ocean_segments+num_land_points],
+               'ocean': range(num_ocean_segments)}
+
+# Build mesh and domain
+domain = create_domain_from_regions(bounding_polygon_gems,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project_GEMS.bounding_maxarea,
+                                    interior_regions=project_GEMS.interior_regions,
+                                    mesh_filename=project_GEMS.meshes,
+                                    use_cache=False, #usually true but changed for testing
+                                    verbose=True)
+print domain.statistics()
+
+domain.starttime(project_GEMS.starttime)
+domain.set_name(project_GEMS.scenario_name)
+domain.set_datadir(project_GEMS.output_run) 
+domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
+
+#-------------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+print 'Setup initial conditions'
+
+# Set the initial stage in the offcoast region only
+if project_GEMS.land_initial_conditions:
+    IC = Polygon_function(project_GEMS.land_initial_conditions,
+                          default=project_GEMS.tide,
+                          geo_reference=domain.geo_reference)
+else:
+    IC = 0
+domain.set_quantity('stage', IC, use_cache=True, verbose=True)
+domain.set_quantity('friction', project_GEMS.friction) 
+domain.set_quantity('elevation', 
+                    filename=project_GEMS.combined_elevation+'.pts',
+                    use_cache=True,
+                    verbose=True,
+                    alpha=project_GEMS.alpha)
+
+#-------------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------------
+
+print 'Set boundary - available tags:', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bt = Transmissive_stage_zero_momentum_boundary(domain)
+Bd = Dirichlet_boundary([project_GEMS.tide, 0, 0])
+Bf = Field_boundary(project_GEMS.event_sts+'.sts',
+                    domain, mean_stage=project_GEMS.tide,
+                    time_thinning=1,
+                    default_boundary=Dirichlet_boundary([0, 0, 0]),
+                    boundary_polygon=bounding_polygon_gems,                    
+                    use_cache=True,
+                    verbose=True)
+
+domain.set_boundary({'back': Br,
+                     'side': Bt,
+                     'ocean': Bf}) 
+
+#-------------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
+t0 = time.time()
+for t in domain.evolve(yieldstep=project_GEMS.yieldstep, 
+                       finaltime=project_GEMS.finaltime,
+                       skip_initial_step=False): 
+    print domain.timestepping_statistics()
+    print domain.boundary_statistics(tags='ocean')
+
+print 'Simulation took %.2f seconds' % (time.time()-t0)
+
Index: /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/setup_model_GEMS.py
===================================================================
--- /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/setup_model_GEMS.py	(revision 7628)
+++ /anuga_work/production/mandurah_storm_surge_2009/gems_comparison/setup_model_GEMS.py	(revision 7628)
@@ -0,0 +1,160 @@
+"""
+Module to import the project_GEMS.py 'configuration' file and perform
+sanity checks plus a quick check on mesh generation.
+
+Also callable as a stand-alone program, mainly to view the results
+of the mesh generation.
+"""
+
+import os
+
+from os.path import join, exists
+from anuga.utilities.polygon import read_polygon, number_mesh_triangles
+
+import project_GEMS
+
+
+#-------------------------------------------------------------------------------
+# Sanity checks - check that required files/directories, etc, exist.
+#-------------------------------------------------------------------------------
+
+# flag - we check many things and then don't proceed if anything wrong
+sanity_error = False               # checked at bottom of this file
+
+# Test that environment variables are defined.
+if os.getenv(project_GEMS.ENV_INUNDATIONHOME) is None:
+	print "Environment variable '%s' is not set" % project_GEMS.ENV_INUNDATIONHOME
+	sanity_error = True
+
+
+#-------------------------------------------------------------------------------
+# Check filename strings that MUST be set
+#-------------------------------------------------------------------------------
+
+if not project_GEMS.gems_order_filename:
+	print ("Sorry, ordering file parameter 'gems_order_filename' "
+		   "isn't set")
+	sanity_error = True
+
+if not project_GEMS.landward_boundary_filename:
+	print ("Sorry, landward bounding points parameter "
+		   "'landward_boundary_filename' isn't set")
+	sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Directory Structure
+#-------------------------------------------------------------------------------
+
+# check folders generated from environment variables.
+if not exists(project_GEMS.home):
+	print "Sorry, data directory '%s' doesn't exist" % project_GEMS.home
+	sanity_error = True
+	
+if not exists(project_GEMS.anuga_folder):
+	print "Sorry, ANUGA directory '%s' doesn't exist" % project_GEMS.anuga_folder
+	sanity_error = True
+
+if not exists(project_GEMS.topographies_folder):
+	print ("Sorry, topo directory '%s' doesn't exist"
+		   % project_GEMS.topographies_folder)
+	sanity_error = True
+
+if not exists(project_GEMS.polygons_folder):
+	print ("Sorry, polygon directory '%s' doesn't exist"
+		   % project_GEMS.polygons_folder)
+	sanity_error = True
+
+if not exists(project_GEMS.boundaries_folder):
+	print ("Sorry, boundaries directory '%s' doesn't exist"
+		   % project_GEMS.boundaries_folder)
+	sanity_error = True
+
+if not exists(project_GEMS.output_folder):
+	print "Sorry, outputs directory '%s' doesn't exist" % project_GEMS.output_folder
+	sanity_error = True
+
+if not exists(project_GEMS.gauges_folder):
+	print "Sorry, gauges directory '%s' doesn't exist" % project_GEMS.gauges_folder
+	sanity_error = True
+
+if not exists(project_GEMS.event_folder):
+	print ("Sorry, you must generate event %s."
+		   % project_GEMS.event)
+	sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Determine type of run
+#-------------------------------------------------------------------------------
+
+if project_GEMS.setup == 'trial':
+	print 'trial'
+	project_GEMS.scale_factor = 100
+	project_GEMS.time_thinning = 96
+	project_GEMS.yieldstep = 240
+elif project_GEMS.setup == 'basic': 
+	print 'basic'
+	project_GEMS.scale_factor = 4
+	project_GEMS.time_thinning = 12
+	project_GEMS.yieldstep = 120
+elif project_GEMS.setup == 'final': 
+	print 'final'
+	project_GEMS.scale_factor = 1
+	project_GEMS.time_thinning = 4
+	project_GEMS.yieldstep = 60
+elif project_GEMS.setup == 'storm_surge_final':
+	print 'storm surge final'
+	project_GEMS.scale_factor = 1
+	project_GEMS.time_thinning = 4
+	project_GEMS.yieldstep = 720	
+else:
+	print ("Sorry, you must set the 'setup' variable to one of:"
+		   '   trial - coarsest mesh, fast\n'
+		   '   basic - coarse mesh\n'
+		   '   final - fine mesh, slowest\n'
+		   '\n'
+		   "'setup' was set to '%s'" % project_GEMS.setup)
+	sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Check for errors detected above.
+#-------------------------------------------------------------------------------
+
+if sanity_error:
+	msg = 'You must fix the above errors before continuing.'
+	raise Exception, msg
+	
+#-------------------------------------------------------------------------------
+# Reading elevation clip polygon
+#-------------------------------------------------------------------------------
+
+project_GEMS.elevation_clip_box = read_polygon(join(project_GEMS.topographies_folder,
+											 project_GEMS.elevation_clip_box_filename))
+
+#-------------------------------------------------------------------------------
+# Reading polygons and creating interior regions
+#-------------------------------------------------------------------------------
+
+# Create list of land polygons with initial conditions
+project_GEMS.land_initial_conditions = []
+for filename, MSL in project_GEMS.land_initial_conditions_filename:
+	polygon = read_polygon(join(project_GEMS.polygons_folder, filename))
+	project_GEMS.land_initial_conditions.append([polygon, MSL])
+
+# Create list of interior polygons with scaling factor
+project_GEMS.interior_regions = []
+for filename, maxarea in project_GEMS.interior_regions_data:
+	polygon = read_polygon(join(project_GEMS.polygons_folder, filename))
+	project_GEMS.interior_regions.append([polygon,
+									 maxarea*project_GEMS.scale_factor])
+
+# Initial bounding polygon for data clipping 
+project_GEMS.bounding_polygon = read_polygon(join(project_GEMS.polygons_folder,
+											 project_GEMS.bounding_polygon_filename))
+project_GEMS.bounding_maxarea = project_GEMS.bounding_polygon_maxarea*project_GEMS.scale_factor
+
+# Estimate the number of triangles                     
+trigs_min = number_mesh_triangles(project_GEMS.interior_regions,
+								  project_GEMS.bounding_polygon,
+								  project_GEMS.bounding_maxarea)
+
+print 'min estimated number of triangles', trigs_min
