Index: /anuga_work/production/Broome_2009/build_elevation.py
===================================================================
--- /anuga_work/production/Broome_2009/build_elevation.py	(revision 7179)
+++ /anuga_work/production/Broome_2009/build_elevation.py	(revision 7179)
@@ -0,0 +1,100 @@
+"""Build the elevation data to run a tsunami inundation scenario 
+for Busselton, WA, 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 import project   # Definition of file names and polygons
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+copy_code_files(project.output_build,__file__, 
+                os.path.dirname(project.__file__)+os.sep+\
+                project.__name__+'.py' )
+start_screen_catcher(project.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.bounding_polygon
+print 'project.combined_elevation_basename', project.combined_elevation_basename
+
+# Create Geospatial data from ASCII files
+geospatial_data = {}
+if not project.ascii_grid_filenames == []:
+    for filename in project.ascii_grid_filenames:
+        absolute_filename = join(project.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.bounding_polygon)
+
+# Create Geospatial data from TXT files
+if not project.point_filenames == []:
+    for filename in project.point_filenames:
+        absolute_filename = join(project.topographies_folder, filename)
+        G_points = Geospatial_data(file_name=absolute_filename,
+                                                    verbose=True)
+        print 'Clip geospatial object'
+        geospatial_data[filename] = G_points.clip(project.bounding_polygon)
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+
+print 'Add geospatial objects' # except', project.offshore_name5
+G = None
+for key in geospatial_data:
+    #if key == project.point_filenames[4] or key == project.ascii_filenames[1]:
+    #    # Skip these files for now
+    #    continue
+    
+    G += geospatial_data[key]
+        
+#print 'Clip combined geospatial data'
+##G_clip = G.clip_outside(project.poly_aoi1)
+##G_all = G_clip + geospatial_data[project.point_filenames[4]] 
+#G_clipped = G_all.clip(project.poly_all)
+#G_clip = G.clip(project.bounding_polygon)
+
+
+print 'Export combined DEM file'
+G.export_points_file(project.combined_elevation + '.pts')
+print 'Do txt version too'
+# Use for comparision in ARC
+G.export_points_file(project.combined_elevation + '.txt')
+
Index: /anuga_work/production/Broome_2009/build_urs_boundary.py
===================================================================
--- /anuga_work/production/Broome_2009/build_urs_boundary.py	(revision 7179)
+++ /anuga_work/production/Broome_2009/build_urs_boundary.py	(revision 7179)
@@ -0,0 +1,167 @@
+"""
+Replacement function that wraps the legacy function urs2sts()
+and passes in the data from the EventSelection <event>.list file.
+"""
+
+import os
+import os.path
+from time import localtime, strftime, gmtime
+
+import project
+
+from Scientific.IO.NetCDF import NetCDFFile
+import Numeric as num
+
+from anuga.shallow_water.data_manager import urs2sts
+
+import Numeric as num
+
+
+#-------------------------------------------------------------------------------
+# Get gauges (timeseries of index points)
+#-------------------------------------------------------------------------------
+def get_sts_gauge_data(filename, verbose=False):
+    print 'get_sts_gauge_data: filename=%s' % filename
+    fid = NetCDFFile(filename+'.sts', 'r')      #Open existing file for read
+    permutation = fid.variables['permutation'][:]
+    x = fid.variables['x'][:] + fid.xllcorner   #x-coordinates of vertices
+    y = fid.variables['y'][:] + fid.yllcorner   #y-coordinates of vertices
+    points = num.transpose(num.asarray([x.tolist(), y.tolist()]))
+    time = fid.variables['time'][:] + fid.starttime
+    elevation = fid.variables['elevation'][:]
+       
+    basename = 'sts_gauge'
+    quantity_names = ['stage', 'xmomentum', 'ymomentum']
+    quantities = {}
+    for i, name in enumerate(quantity_names):
+        quantities[name] = fid.variables[name][:]
+
+    #---------------------------------------------------------------------------
+    # Get maximum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+     
+    maxname = 'max_sts_stage.csv'
+    print 'get_sts_gauge_data: maxname=%s' % maxname
+    fid_max = open(os.path.join(project.event_folder, maxname), 'w')
+    fid_max.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage)))
+      
+    #---------------------------------------------------------------------------
+    # Get minimum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+
+    minname = 'min_sts_stage.csv'
+    fid_min = open(os.path.join(project.event_folder, minname), 'w')
+    fid_min.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)))
+
+        out_file = os.path.join(project.event_folder,
+                                basename+'_'+str(index)+'.csv')
+        fid_sts = open(out_file, 'w')
+        fid_sts.write('time, stage, xmomentum, ymomentum \n')
+
+        #-----------------------------------------------------------------------
+        # End of the get gauges
+        #-----------------------------------------------------------------------
+        for k in range(len(time)-1):
+            fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
+                          % (time[k], stage[k], xmomentum[k], ymomentum[k]))
+
+        fid_sts.close()      
+
+    fid.close()
+
+    return quantities,elevation,time
+
+
+##
+# @brief Build boundary STS files from one or more MUX files.
+# @param event_file Name of mux meta-file or single mux stem.
+# @param output_dir Directory to write STS data to.
+# @note 'event_file' is produced by EventSelection.
+def build_urs_boundary(event_file, output_dir):
+    '''Build a boundary STS file from a set of MUX files.'''
+
+    # if we are using an EventSelection multi-mux file
+    if project.multi_mux:
+        # get the mux+weight data from the meta-file (in <boundaries>)
+        mux_event_file = os.path.join(project.event_folder, event_file)
+        print 'using multi-mux file', mux_event_file
+        try:
+            fd = open(mux_event_file, 'r')
+            mux_data = fd.readlines()
+            fd.close()
+        except IOError, e:
+            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
+            raise Exception, msg
+        except:
+            raise
+
+        # first line of file is # filenames+weight in rest of file
+        num_lines = int(mux_data[0].strip())
+        mux_data = mux_data[1:]
+        print 'number of sources %d' % num_lines
+
+        # quick sanity check on input mux meta-file
+        if num_lines != len(mux_data):
+            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
+                   % (event_file, len(mux_data), num_lines))
+            raise Exception, msg
+
+        # Create filename and weights lists.
+        # Must chop GRD filename just after '*.grd'.
+        mux_filenames = []
+        for line in mux_data:
+            muxname = line.strip().split()[0]
+            split_index = muxname.index('.grd')
+            muxname = muxname[:split_index+len('.grd')]
+            muxname = os.path.join(project.mux_data_folder, muxname)
+            mux_filenames.append(muxname)
+        
+        mux_weights = [float(line.strip().split()[1]) for line in mux_data]
+
+        # Call legacy function to create STS file.
+        print 'creating sts file'
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=project.urs_order,
+                weights=mux_weights,
+                zone=project.zone,
+                verbose=True)
+    else:                           # a single mux stem file, assume 1.0 weight
+        mux_file = os.path.join(project.event_folder, event_file)
+        mux_filenames = [mux_file]
+        print 'using single-mux file', mux_file
+
+        weight_factor = 1.0
+        mux_weights = weight_factor*num.ones(len(mux_filenames), num.Float)
+            
+        order_filename = project.urs_order
+
+        print 'reading', order_filename
+        # Create ordered sts file
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=order_filename,
+                weights=mux_weights,
+                verbose=True)
+
+    # report on progress so far
+    sts_file = os.path.join(project.event_folder, project.scenario_name)
+    quantities, elevation, time = get_sts_gauge_data(sts_file, verbose=False)
+    print len(elevation), len(quantities['stage'][0,:])
+
+    
+
Index: /anuga_work/production/Broome_2009/file_length.py
===================================================================
--- /anuga_work/production/Broome_2009/file_length.py	(revision 7179)
+++ /anuga_work/production/Broome_2009/file_length.py	(revision 7179)
@@ -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/Broome_2009/project.py
===================================================================
--- /anuga_work/production/Broome_2009/project.py	(revision 7179)
+++ /anuga_work/production/Broome_2009/project.py	(revision 7179)
@@ -0,0 +1,254 @@
+"""
+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 = 'broome'
+scenario_folder = 'broome_tsunami_scenario_2009'
+
+#-------------------------------------------------------------------------------
+# Initial Conditions
+#-------------------------------------------------------------------------------
+
+# Model specific parameters.
+# One or all can be changed each time the run_model script is executed
+tide = 0.0                # difference between MSL and HAT in metres
+zone = 51               # specify UTM zone of model
+event_number = 27279    # 58280, 64477 the event number or the mux file name
+alpha = 0.1             # smoothing parameter for mesh
+friction=0.01           # manning's friction coefficient
+starttime=0             # start time for simulation
+finaltime=1000         # final time for simulation
+
+setup = 'trial'         # This can be one of three values
+                        #    trial - coarsest mesh, fast
+                        #    basic - coarse mesh
+                        #    final - fine mesh, slowest
+
+#-------------------------------------------------------------------------------
+# 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_number]
+
+#-------------------------------------------------------------------------------
+# Input Data
+#-------------------------------------------------------------------------------
+
+# ELEVATION DATA
+# Used in build_elevation.py
+ascii_grid_filenames = ['town_topo_10m',
+                        'cable_250m',
+                        'inferred_north',
+                        'inferred_south',
+                        'north_250m',
+                        'other_topo_250m',
+                        'south_250m']
+point_filenames = ['Broome_coastline.txt','Broome_Bathymetry.txt']
+
+### Add csv header list to all files in point_filenames
+##headerlist = ['x', 'y', 'elevation'] 
+##for f in point_filenames:
+##    add_csv_header(join(topographies_folder, f), headerlist)
+
+# BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
+# Used in build_elevation.py
+# Format for points easting,northing (no header)
+bounding_polygon_filename = 'bounding_polygon.csv'
+bounding_polygon_maxarea = 1000000
+
+# INTERIOR REGIONS -  for designing the mesh
+# Used in run_model.py
+# Format for points easting,northing (no header)                    
+interior_regions_data = [['Shallow_Water_MH.csv', 1500],
+                         ['AoS_MH', 800],
+                         ['AoI_MH', 500],
+                         ['AoS_north_MH', 800],
+                         ['AoI_north_MH', 500],
+                         ['bay_coast_MH', 1500],
+                         ['bay_coast_small_MH', 800]]
+PriorityArea_filename = 'PriorityAreas.csv'
+    
+# 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_mainland.csv'
+
+# GAUGES - for creating timeseries at a specific point
+# Used in get_timeseries.py.  
+# Format easting,northing,name,elevation (with header)
+gauges_filename = ''
+# BUILDINGS EXPOSURE - for identifying inundated houses
+# Used in run_building_inundation.py
+# Format latitude,longitude etc (geographic)
+building_exposure_filename = 'busselton_res_clip.csv' # from NEXIS
+
+# AREA OF IMAGES - Extent of each image to find out highest runup
+# Header - easting,northing,id,value
+# Used in get_runup.py
+images_filename = ''
+
+# BOUNDING POLYGON - used in build_boundary.py and run_model.py respectively
+# NOTE: when files are put together the points must be in sequence
+# For ease go clockwise!
+# Check the run_model.py for boundary_tags
+
+# Thinned ordering file from Hazard Map (geographic)
+# Format is index,latitude,longitude (with header)
+urs_order_filename = 'thinned_boundary_ordering.csv'
+
+# Landward bounding points
+# Format easting,northing (no header)
+landward_boundary_filename = 'landward_boundary_MH.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
+#-------------------------------------------------------------------------------
+
+# ASCII export grid for Bruny
+##xminBruny = 523900
+##xmaxBruny = 533200
+##yminBruny = 5204300
+##ymaxBruny = 5213100
+
+
+################################################################################
+################################################################################
+####         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'
+
+# Path to MUX data
+ENV_MUXHOME = 'MUXHOME'
+
+#-------------------------------------------------------------------------------
+# 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 + '_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
+muxhome = os.getenv(ENV_MUXHOME)
+    
+# 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_number))
+
+# MUX data files
+# Directory containing the MUX data files to be used with EventSelection.
+mux_data_folder = join(muxhome, 'mux')
+
+#-------------------------------------------------------------------------------
+# 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 urs order points, used within build_urs_boundary.py
+if urs_order_filename:
+    urs_order = join(boundaries_folder, urs_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) + '.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)
+
+#Multiple polygons in one CSV file to make internal polygons
+if not PriorityArea_filename == None:
+    PriorityAreas = join(polygons_folder, PriorityArea_filename)
+
+
Index: /anuga_work/production/Broome_2009/run_model.py
===================================================================
--- /anuga_work/production/Broome_2009/run_model.py	(revision 7179)
+++ /anuga_work/production/Broome_2009/run_model.py	(revision 7179)
@@ -0,0 +1,164 @@
+"""Run a tsunami inundation scenario for Busselton, WA, Australia.
+
+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 tsunami is generated through an sts file from build_boundary.py.
+
+Input: sts file (build_boundary.py 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
+"""
+
+#------------------------------------------------------------------------------
+# 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 Numeric 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 import project
+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.output_run, __file__, 
+                os.path.join(os.path.dirname(project.__file__),
+                             project.__name__+'.py'))
+start_screen_catcher(project.output_run, 0, 1)
+
+#-------------------------------------------------------------------------------
+# Create the computational domain 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
+#-------------------------------------------------------------------------------
+
+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
+event_sts = create_sts_boundary(project.event_sts)
+
+# Reading the landward defined points, this incorporates the original clipping
+# polygon minus the 100m contour
+landward_boundary = read_polygon(project.landward_boundary)
+
+# Combine sts polyline with landward points
+bounding_polygon_sts = event_sts + landward_boundary
+
+# Number of boundary segments
+num_ocean_segments = len(event_sts) - 1
+# Number of landward_boundary points
+num_land_points = file_length(project.landward_boundary)
+
+# Boundary tags refer to project.landward_boundary
+# 4 points equals 5 segments start at N
+boundary_tags={'back': range(num_ocean_segments+1,
+                             num_ocean_segments+num_land_points-1), #last segment will be side
+               'side': [num_ocean_segments,
+                        num_ocean_segments+num_land_points-1, # three side segments total
+                        num_ocean_segments+num_land_points],
+               'ocean': range(num_ocean_segments)}
+
+# Build mesh and domain
+domain = create_domain_from_regions(bounding_polygon_sts,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project.bounding_maxarea,
+                                    interior_regions=project.interior_regions,
+                                    mesh_filename=project.meshes,
+                                    use_cache=True,
+                                    verbose=True)
+print domain.statistics()
+
+domain.set_name(project.scenario_name)
+domain.set_datadir(project.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.land_initial_conditions:
+    IC = Polygon_function(project.land_initial_conditions,
+                          default=project.tide,
+                          geo_reference=domain.geo_reference)
+else:
+    IC = 0
+domain.set_quantity('stage', IC, use_cache=True, verbose=True)
+domain.set_quantity('friction', project.friction) 
+domain.set_quantity('elevation', 
+                    filename=project.combined_elevation+'.pts',
+                    use_cache=True,
+                    verbose=True,
+                    alpha=project.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.tide, 0, 0])
+Bf = Field_boundary(project.event_sts+'.sts',
+                    domain, mean_stage=project.tide,
+                    time_thinning=1,
+                    default_boundary=Dirichlet_boundary([0, 0, 0]),
+                    boundary_polygon=bounding_polygon_sts,                    
+                    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.yieldstep, 
+                       finaltime=project.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/Broome_2009/setup_model.py
===================================================================
--- /anuga_work/production/Broome_2009/setup_model.py	(revision 7179)
+++ /anuga_work/production/Broome_2009/setup_model.py	(revision 7179)
@@ -0,0 +1,198 @@
+"""
+Module to import the project.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
+from anuga.shallow_water.data_manager import csv2polygons
+import project
+
+
+#-------------------------------------------------------------------------------
+# 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.ENV_INUNDATIONHOME) is None:
+    print "Environment variable '%s' is not set" % project.ENV_INUNDATIONHOME
+    sanity_error = True
+
+if os.getenv(project.ENV_MUXHOME) is None:
+    print "Environment variable '%s' is not set" % project.ENV_MUXHOME
+    sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Check filename strings that MUST be set
+#-------------------------------------------------------------------------------
+
+if not project.urs_order_filename:
+    print ("Sorry, thinned ordering file parameter 'urs_order_filename' "
+           "isn't set")
+    sanity_error = True
+
+if not project.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.home):
+    print "Sorry, data directory '%s' doesn't exist" % project.home
+    sanity_error = True
+    
+if not exists(project.muxhome):
+    print "Sorry, MUX directory '%s' doesn't exist" % project.muxhome
+    sanity_error = True
+    
+if not exists(project.anuga_folder):
+    print "Sorry, ANUGA directory '%s' doesn't exist" % project.anuga_folder
+    sanity_error = True
+
+if not exists(project.topographies_folder):
+    print ("Sorry, topo directory '%s' doesn't exist"
+           % project.topographies_folder)
+    sanity_error = True
+
+if not exists(project.polygons_folder):
+    print ("Sorry, polygon directory '%s' doesn't exist"
+           % project.polygons_folder)
+    sanity_error = True
+
+if not exists(project.boundaries_folder):
+    print ("Sorry, boundaries directory '%s' doesn't exist"
+           % project.boundaries_folder)
+    sanity_error = True
+
+if not exists(project.output_folder):
+    print "Sorry, outputs directory '%s' doesn't exist" % project.output_folder
+    sanity_error = True
+
+if not exists(project.gauges_folder):
+    print "Sorry, gauges directory '%s' doesn't exist" % project.gauges_folder
+    sanity_error = True
+
+if not exists(project.mux_data_folder):
+    print "Sorry, mux data directory '%s' doesn't exist" % project.mux_data_folder
+    sanity_error = True
+
+# if multi_mux is True, check if multi-mux file exists
+if project.multi_mux:
+    if not exists(project.mux_input):
+        print ("Sorry, MUX input file '%s' doesn't exist"
+               % project.mux_input)
+        sanity_error = True
+
+#-----
+# If this directory don't exist, EventSelection hasn't been run.
+#-----
+
+if not exists(project.event_folder):
+    print ("Sorry, you must generate event %s with EventSelection."
+           % project.event_number)
+    sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Determine type of run
+#-------------------------------------------------------------------------------
+
+if project.setup == 'trial':
+    print 'trial'
+    project.scale_factor = 100
+    project.time_thinning = 96
+    project.yieldstep = 240
+elif project.setup == 'basic': 
+    print 'basic'
+    project.scale_factor = 4
+    project.time_thinning = 12
+    project.yieldstep = 120
+elif project.setup == 'final': 
+    print 'final'
+    project.scale_factor = 1
+    project.time_thinning = 4
+    project.yieldstep = 60
+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.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 polygons and creating interior regions
+#-------------------------------------------------------------------------------
+
+# Create list of land polygons with initial conditions
+project.land_initial_conditions = []
+# if it's a list, then it's a list of land condition filenames
+# else it's a string - a single file, multiple land conditions
+if isinstance(project.land_initial_conditions_filename, list):
+    for filename, MSL in project.land_initial_conditions_filename:
+        polygon = read_polygon(join(project.polygons_folder, filename))
+        project.land_initial_conditions.append([polygon, MSL])
+elif isinstance(project.land_initial_conditions_filename, basestring):
+    polygons, MSL = csv2polygons(join(project.polygons_folder,
+                                      project.land_initial_conditions_filename))
+    for i, key in enumerate(polygons):
+        if i%100==0: print i
+        poly = polygons[key]
+        land = float(MSL[key])
+        project.land_initial_conditions.append([poly, land])
+else:
+    msg = ('project.land_initial_conditions_filename must be a list or '
+           'string, got %s'
+           % type(project.land_initial_conditions_filename))
+    raise Exception, msg
+# Create list of interior polygons with scaling factor
+project.interior_regions = []
+
+if project.PriorityArea_filename is not None:
+    polygons, maxareas = csv2polygons(project.PriorityAreas)
+    print 'Creating %d internal polygons' % len(polygons)
+    #def create_polygon_function(polygons, geo_reference=None):
+    project.interior_regions = []
+    for i, key in enumerate(polygons):
+        if i%100==0: print i
+        poly = polygons[key]
+        maxarea = float(maxareas[key])
+        project.interior_regions.append([poly,
+                                         maxarea*project.scale_factor])
+
+for filename, maxarea in project.interior_regions_data:
+    polygon = read_polygon(join(project.polygons_folder, filename))
+    project.interior_regions.append([polygon,
+                                     maxarea*project.scale_factor])
+
+# Initial bounding polygon for data clipping 
+project.bounding_polygon = read_polygon(join(project.polygons_folder,
+                                             project.bounding_polygon_filename))
+project.bounding_maxarea = project.bounding_polygon_maxarea*project.scale_factor
+
+# Estimate the number of triangles                     
+trigs_min = number_mesh_triangles(project.interior_regions,
+                                  project.bounding_polygon,
+                                  project.bounding_maxarea)
+
+print 'min estimated number of triangles', trigs_min
