Index: /anuga_core/documentation/user_manual/demos/cairns/ExportResults.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/cairns/ExportResults.py	(revision 3979)
+++ /anuga_core/documentation/user_manual/demos/cairns/ExportResults.py	(revision 3979)
@@ -0,0 +1,48 @@
+import project, os
+import sys
+
+from anuga.shallow_water.data_manager import sww2dem
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
+from os import sep
+
+time_dir = '20061101_214314' 
+
+directory = project.outputdir
+name = directory + time_dir +sep + 'source'
+
+print 'output dir:', name
+which_var = 2
+if which_var == 0:  # Stage
+    outname = name + '_stage'
+    quantityname = 'stage'
+
+if which_var == 1:  # Absolute Momentum
+    outname = name + '_momentum'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  #Absolute momentum
+
+if which_var == 2:  # Depth
+    outname = name + '_depth_bruny'
+    quantityname = 'stage-elevation'  #Depth
+
+if which_var == 3:  # Speed
+    outname = name + '_speed'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-30)'  #Speed
+
+if which_var == 4:  # Elevation
+    outname = name + '_elevation'
+    quantityname = 'elevation'  #Elevation
+
+print 'start sww2dem'
+
+sww2dem(name, basename_out = outname,
+            quantity = quantityname,
+            cellsize = 25,      # would prefer this at 25
+            # define region for viz purposes
+            easting_min = project.eastingmin,
+            easting_max = project.eastingmax,
+            northing_min = project.northingmin,
+            northing_max = project.northingmax,        
+            reduction = max, #this is because we want max quantityname
+            verbose = True,
+            format = 'asc')
+
Index: /anuga_core/documentation/user_manual/demos/cairns/GetTimeseries.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/cairns/GetTimeseries.py	(revision 3979)
+++ /anuga_core/documentation/user_manual/demos/cairns/GetTimeseries.py	(revision 3979)
@@ -0,0 +1,43 @@
+"""
+Generate time series of nominated "gauges" 
+Note, this script will only work if pylab is installed on the platform
+
+Inputs:
+
+production dirs: dictionary of production directories with a
+                 association to that simulation run, eg high tide,
+                 magnitude, etc.
+                    
+Outputs:
+
+* figures stored in same directory as sww file
+* time series data stored in csv files in same directory as sww file
+* elevation at nominated gauges (elev_output)
+"""
+
+from os import getcwd, sep, altsep, mkdir, access, F_OK
+import project
+from anuga.abstract_2d_finite_volumes.util import sww2timeseries
+
+# nominate directory location of sww file with associated attribute
+production_dirs = {'slump': 'Slump',
+                   'fixedwave': 'Fixed Wave'}
+
+# Generate figures
+swwfiles = {}
+for label_id in production_dirs.keys():
+    file_loc = project.outputdir + label_id + sep
+    swwfile = file_loc + project.basename + production_dirs[label_id] + '.sww'
+    swwfiles[swwfile] = label_id
+
+texname, elev_output = sww2timeseries(swwfiles,
+                                      project.gauge_filename,
+                                      production_dirs,
+                                      report = False,
+                                      reportname = '',
+                                      plot_quantity = ['stage', 'speed'],
+                                      surface = False,
+                                      time_min = None,
+                                      time_max = None,
+                                      title_on = True,
+                                      verbose = True)
Index: /anuga_core/documentation/user_manual/demos/cairns/project.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/cairns/project.py	(revision 3979)
+++ /anuga_core/documentation/user_manual/demos/cairns/project.py	(revision 3979)
@@ -0,0 +1,94 @@
+# -*- 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, \
+                                    polygon_area, is_inside_polygon
+
+if sys.platform == 'win32':
+    home = getenv('INUNDATIONHOME')
+    user = getenv('USERPROFILE')
+
+else:    
+    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+\
+                  '2'+sep+'ramp'+sep+'risk_assessment_methods_project'+\
+                  sep+'inundation')     
+    user = getenv('LOGNAME')
+    print 'USER:', user
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
+
+#Making assumptions about the location of scenario data
+state = 'queensland'
+scenario_dir_name = 'cairns_tsunami_scenario_demo'
+
+# onshore data provided by WA DLI
+demname = 'cairns'
+
+#swollen/ all data output
+basename = 'source'
+codename = 'project.py'
+
+#Derive subdirectories and filenames
+local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
+meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
+datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
+gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
+outputtimedir = outputdir + local_time + sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+
+gauge_filename = gaugedir + 'gauges.csv'
+codedir = getcwd()+sep                            
+codedirname = codedir + 'project.py'
+#meshname = outputtimedir + 'mesh_' + basename
+
+# Necessary if using point datasets, rather than grid
+dem_name = datadir + demname
+
+###############################
+# Domain definitions
+###############################
+
+# bounding polygon for study area
+polyAll = read_polygon(polygondir+'extent_v2.csv')
+
+print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 356300
+eastingmax = 398500
+northingmin = 8107850
+northingmax = 8148700
+
+###############################
+# Interior region definitions
+###############################
+
+# broome digitized polygons
+poly_cairns = read_polygon(polygondir+'cairns.csv')
+poly_island0 = read_polygon(polygondir+'islands.csv')
+poly_island1 = read_polygon(polygondir+'islands1.csv')
+poly_island2 = read_polygon(polygondir+'islands2.csv')
+poly_island3 = read_polygon(polygondir+'islands3.csv')
+poly_shallow = read_polygon(polygondir+'shallow.csv')
+
+plot_polygons([polyAll,poly_cairns,poly_island0,poly_island1,\
+               poly_island2,poly_island3,poly_shallow],\
+              'boundingpoly',verbose=False)
+
+#slump_origin = [633361, 8129700]
+#slump_depth = 986
+slump_origin = [451871, 8128376] # move onto the continental shelf, depth = 500
+slump_depth = 500.
Index: /anuga_core/documentation/user_manual/demos/cairns/runcairns.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/cairns/runcairns.py	(revision 3979)
+++ /anuga_core/documentation/user_manual/demos/cairns/runcairns.py	(revision 3979)
@@ -0,0 +1,223 @@
+"""Script for running a tsunami inundation scenario for Broome, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a tsunami wave generated by MOST.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
+"""
+
+#-------------------------------------------------------------------------------
+# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+from shutil import copy
+from os import mkdir, access, F_OK
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary, Time_boundary, File_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+#-------------------------------------------------------------------------------
+# Define scenario as either slump or fixed_wave.
+#-------------------------------------------------------------------------------
+scenario = 'slump' # 'fixedwave'
+mkdir (scenario)
+basename = scenario + project.basename
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+
+# creates copy of code in output dir if dir doesn't exist
+if access(project.outputtimedir,F_OK) == 0 :
+    mkdir (project.outputtimedir)
+copy (project.codedirname, project.outputtimedir + project.codename)
+copy (project.codedir + 'runcairns.py', project.outputtimedir + 'runcairns2.py')
+print'output dir', project.outputtimedir
+
+#normal screen output is stored in 
+screen_output_name = project.outputtimedir + "screen_output.txt"
+screen_error_name = project.outputtimedir + "screen_error.txt"
+
+#used to catch screen output to file
+sys.stdout = Screen_Catcher(screen_output_name)
+sys.stderr = Screen_Catcher(screen_error_name)
+
+print 'USER:    ', project.user
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+#-------------------------------------------------------------------------------
+
+# filenames
+dem_name = project.dem_name
+meshname = project.meshname+'.msh'
+
+# creates DEM from asc data
+convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+dem2pts(dem_name, use_cache=True, verbose=True)
+
+#----------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with a tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+remainder_res = 10000000
+islands_res = 100000
+cairns_res = 100000
+shallow_res = 500000
+interior_regions = [[project.poly_cairns, cairns_res],
+                    [project.poly_island0, islands_res],
+                    [project.poly_island1, islands_res],
+                    [project.poly_island2, islands_res],
+                    [project.poly_island3, islands_res],
+                    [project.poly_shallow, shallow_res]]
+
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project.polyAll,
+           {'boundary_tags': {'top': [0], 'ocean_east': [1],
+                              'bottom': [2], 'onshore': [3]},
+           'maximum_triangle_area': remainder_res,
+           'filename': meshname,
+           'interior_regions': interior_regions},
+          verbose = True, evaluate=False)
+print 'created mesh'
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+domain = Domain(meshname, use_cache = True, verbose = True)
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+print domain.statistics()
+
+domain.set_name(basename) #domain.set_name(project.basename)
+domain.set_datadir(scenario) #domain.set_datadir(project.outputtimedir)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
+
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+tide = 0.0
+domain.set_quantity('stage', tide)
+domain.set_quantity('friction', 0.0) 
+domain.set_quantity('elevation', 
+                    filename = project.dem_name + '.pts',
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.1
+                    )
+
+#-------------------------------------------------------------------------------                                 
+# Setup information for slump scenario (to be applied 1 min into simulation
+#-------------------------------------------------------------------------------
+if scenario == 'slump':
+    from anuga.shallow_water.smf import slump_tsunami  # Function for submarine slump
+    tsunami_source = slump_tsunami(length=15000.0,
+                                   depth=project.slump_depth,
+                                   slope=6.0,
+                                   thickness=250.0, 
+                                   radius=3330,
+                                   dphi=0.23,
+                                   x0=project.slump_origin[0], 
+                                   y0=project.slump_origin[1], 
+                                   alpha=0.0, 
+                                   domain=domain,
+                                   verbose=True)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+
+# 7 min square wave starting at 1 min, 10m high
+if scenario == 'fixed_wave':
+    Bw = Time_boundary(domain = domain,
+                       f=lambda t: [(60<t<480)*10, 0, 0])
+    domain.set_boundary( {'ocean_east': Bw, 'bottom': Bd, 'onshore': Bd,
+                           'top': Bd} )
+
+# boundary conditions for slump scenario
+if scenario == 'slump':
+    domain.set_boundary( {'ocean_east': Bd, 'bottom': Bd, 'onshore': Bd,
+                          'top': Bd} )
+
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+import time
+t0 = time.time()
+
+from Numeric import allclose
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+if scenario == 'slump':
+    
+    for t in domain.evolve(yieldstep = 30, finaltime = 60): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean_east')     
+        
+        # add slump
+        thisstagestep = domain.get_quantity('stage')
+        if allclose(t, 60):
+            slump = Quantity(domain)
+            slump.set_values(tsunami_source)
+            domain.set_quantity('stage', slump + thisstagestep)
+
+        # save every two mins leading up to wave approaching land
+        for t in domain.evolve(yieldstep = 120, finaltime = 1000, 
+                               skip_initial_step = True):
+            domain.write_time()
+            domain.write_boundary_statistics(tags = 'ocean_east')
+
+        # save every 30 secs as wave starts inundating ashore
+        for t in domain.evolve(yieldstep = 30, finaltime = 10000, 
+                               skip_initial_step = True):
+            domain.write_time()
+            domain.write_boundary_statistics(tags = 'ocean_east')
+
+if scenario == 'fixed_wave':
+
+    # save every two mins leading up to wave approaching land
+    for t in domain.evolve(yieldstep = 120, finaltime = 1000): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean_east')     
+
+        # save every 30 secs as wave starts inundating ashore
+        for t in domain.evolve(yieldstep = 30, finaltime = 10000, 
+                               skip_initial_step = True):
+            domain.write_time()
+            domain.write_boundary_statistics(tags = 'ocean_east')
+            
+print 'That took %.2f seconds' %(time.time()-t0)
