Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/make_anugaviewer_movie.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/make_anugaviewer_movie.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/make_anugaviewer_movie.py	(revision 9336)
@@ -0,0 +1,115 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+
+Script to automate steps in making an anuga-viewer based movie.
+
+This assumes you have made a directory of jpg images with anuga-viewer.
+
+It also assumes the 'convert' program is installed [part of ImageMajik].
+
+#######################
+To produce the images:
+
+1) Start anuga-viewer:
+
+> anuga_viewer mysww.sww
+
+and get the view oriented how you would like for the movie.
+
+NOTE: To cull very shallow depths, you can do things like
+
+> anuga_viewer -hmin 0.00001 mysww.sww
+
+You need to experiment to get a good value of hmin [it is not the physical
+depth]
+
+2) Press 1 to start recording, space to start the animation, then press 1 to
+finish recording, and 3 to save the movie to 'movie.swm'
+
+3) Close anuga-viewer
+
+4) Make a directory with images of the movie, with:
+
+> anuga-viewer -movie MYMOVIEDIR movie.swm
+
+Let it run until finished.
+
+5) Look at the images. You might want to delete the first one (seems to have a
+blue background irrespective of the actual value of trace set in anuga-viewer).
+Don't worry about repeats -- this script will remove them.
+
+Finally you can run the current script, setting the input parameters as
+appropriate
+
+> python makeAnugaViewerMovie.py
+
+"""
+import os
+import sys
+import glob
+import numpy
+import shutil
+
+# INPUT parameters
+
+movie_dir = 'MYMOVIEDIR'  # Directory where jpg images to make movie are stored
+delay = 20  # Delay between movie frames, in (1/100)ths of a second
+
+# You probably don't want to change the following parameters
+
+file_wildcard = 'frame_*.jpg'
+animation_name = 'animation.gif'  # Must end in .gif
+
+# END INPUT parameters
+
+#################################################################
+
+os.chdir(movie_dir)
+
+# File names
+
+alljpg = glob.glob(file_wildcard)
+
+# File name integer
+
+m = len(alljpg[0].split('_'))
+alljpg_int = [i.split('_')[m - 1] for i in alljpg]
+alljpg_int = [int(i.split('.')[0]) for i in alljpg_int]
+
+# File sizes
+
+alljpg_size = [os.stat(i)[6] for i in alljpg]
+
+# Sort file names and sizes
+
+new_inds = numpy.array(alljpg_int).argsort()
+alljpg = [alljpg[new_inds[i]] for i in range(len(new_inds))]
+alljpg_size = [alljpg_size[new_inds[i]] for i in range(len(new_inds))]
+
+# Make a temporary directory
+
+tmp_dir = 'TMP/'
+os.mkdir(tmp_dir)
+
+# Copy all files with unique size to the temporary directory
+# To ensure the order is interpreted correctly, we number the file names
+
+for i in range(len(alljpg)):
+    if i == 0:
+        shutil.copy2(alljpg[i], tmp_dir + str(10000 + i) + '_'
+                     + os.path.basename(alljpg[i]))
+    else:
+        if alljpg_size[i] != alljpg_size[i - 1]:
+            shutil.copy2(alljpg[i], tmp_dir + str(10000 + i) + '_'
+                         + os.path.basename(alljpg[i]))
+
+# Run movie command
+movie_command = 'convert -delay ' + str(delay) + ' TMP/*' \
+    + file_wildcard + ' ' + animation_name
+print movie_command
+os.system(movie_command)
+
+# Clean up directory
+
+shutil.rmtree(tmp_dir)
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/project.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/project.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/project.py	(revision 9336)
@@ -0,0 +1,417 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+
+THE USER MUST EDIT THIS FILE TO DEFINE THE MESH / RESOLUTION / FILE LOCATIONS
+ETC
+
+Setup base data for ANUGA run
+
+"""
+
+
+import glob
+import os
+import time
+import shutil
+import hashlib
+import json
+import numpy
+from anuga.utilities import spatialInputUtil as su
+import anuga
+from anuga_parallel import myid, barrier, send, receive, numprocs
+
+# We need the boundary tags info to correctly cache the mesh
+# from setup_boundary_conditions import boundary_info
+
+import user_functions
+
+# @@@@@@@@@@@@@@@@@#
+# START EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+#####################################################################
+#
+# SCENARIO SETUP
+#
+#####################################################################
+
+# Name for scenario
+#
+
+scenario = 'tsunami1_template'
+
+# Total model runtime in seconds (starts at time=0.)
+#
+
+finaltime = 0.05 * 3600.0
+
+# Model time between output file timesteps in seconds
+#
+
+yieldstep = 20.0
+
+# Spatial Information
+#
+# Give UTM Zone as an integer (negative for Southern Hemisphere)
+# or a proj4string defining the coordinate system (which is more flexible)
+
+projection_information = 47
+
+# projection_information='+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs'
+
+#####################################################################
+#
+# INPUT DATA DEFINING MESH
+#
+#####################################################################
+
+# Bounding polygon with default mesh resolution (as a max triangle area in m^2)
+#
+# Filename MUST be a .shp file, and contain a set of lines making up the
+# boundary polygon The lines each have an attribute giving the boundary tag, in
+# the attribute table column named (the value of boundary_tags_attribute_name)
+#
+# The starting/ending points on each line MUST EXACTLY AGREE so the whole
+# thing makes a polygon.
+#
+# Make it in GIS by first making the bounding polygon, then splitting it at
+# the appropriate locations, then adding the boundary tag attribute.
+#
+# bounding_polygon_file='../anuga/polygons/bounding_polygon.csv'
+
+bounding_polygon_and_tags_file = \
+    '../anuga/boundaries/Patong_bounding_polygon/Patong_bounding_polygon.shp'
+boundary_tags_attribute_name = 'bndryTag'
+
+default_res = 0.5 * 500 ** 2
+
+# Interior regions with maximum triangle area resolutions (m^2)
+#
+# Filenames can be .shp or .csv
+#
+# If this is an empty list, we can instead use breaklines + ptAreas defined
+# spatially
+#
+
+interior_regions_data = [
+    ['../anuga/polygons/inundation_area.csv', 75],
+    ['../anuga/polygons/aoi.csv', 200],
+    ['../anuga/polygons/aos.csv', 900],
+    ['../anuga/polygons/sw_coast_original.csv', 7000],
+    ['../anuga/polygons/building_main.csv', 20],
+    ['../anuga/polygons/building_main_small.csv', 20],
+    ['../anuga/polygons/building_main_south.csv', 19],
+    ['../anuga/polygons/building_saddle.csv', 20],
+]
+
+# BreakLines.
+#
+# List of filenames (can be .shp and/or .csv)
+#
+# breakline_files=glob.glob('../../GIS/Mesh/BreakLines/*.shp')
+
+breakline_files = []
+
+# Riverwalls
+#
+# Filenames MUST all be .csv with at least x,y,z of riverwall, and and optional
+# header describing the riverwallPar
+#
+# riverwall_csv_files=glob.glob('../../GIS/Mesh/RiverWalls/*.csv')
+
+riverwall_csv_files = []
+
+# If breaklines/riverwalls intersect with < this distance between the
+# intersection and an existing point, then the existing point is moved to the
+# intersection. Otherwise the intersection point is added
+
+break_line_intersect_point_movement_threshold = 10.
+
+# Get mesh resolutions = region_point_areas [alternative to interior_polygons +
+# res] ptAreas='../../GIS/Mesh/RegionResolutions/RegionResolutions.shp'
+
+pt_areas = None
+
+######################################################################
+#
+# INITIAL CONDITIONS (minimal)
+#
+# We prefer to define data for initial / boundary conditions in the files
+# setup_initial_conditions
+#         and
+# setup_boundary_conditions
+#
+# However some variables are needed globally, so they can go here
+#
+######################################################################
+
+# Initial 'flat' water surface
+#
+
+tide = 0.
+
+#####################################################################
+#
+# Numerics
+#
+#####################################################################
+
+# ANUGA's flow algorithm
+#
+
+flow_algorithm = 'DE1'
+
+# Local extrapolation and flux updating
+#
+# Only set this to True if you know what you are doing
+#
+
+use_local_extrapolation_and_flux_updating = False
+
+#########################################################################
+#
+# OUTPUT INFORMATION
+#
+#########################################################################
+
+# Output tif details
+#
+# pixle size (m) for output tifs
+
+output_tif_cellsize = 250.
+
+# Polygon for output tif [either a filename, or if None the
+# bounding_polygon is used]
+
+output_tif_bounding_polygon = None
+
+# Storage of max quantities
+#
+# Stage/depth/speed/speedxDepth maxima are stored separately to the SWW file
+# Need to determine whether they are updated every time step or less frequently
+# (less frequently might save computational time ?? )
+# Can also choose to ignore an initial time period
+
+# 1 = Every time step, 2 = every second timestep, etc
+max_quantity_update_frequency = 1
+# Set this to e.g. 3600.*2 to ignore the first 2 hours
+max_quantity_collection_start_time = 0.
+
+# sww file-size options
+#
+# If store_vertices_uniquely=True, then SWW files are larger and vertices are
+# stored uniquely for every triangle. Otherwise, vertex values are averaged
+# before storage, leading to smaller file sizes but a less accurate
+# representation of what the algorithm is doing
+
+store_vertices_uniquely = False
+
+# Checkpointing
+#
+# Need to restructure to allow checkpointing
+# checkpoint_time=30
+
+# @@@@@@@@@@@@@@@@@#
+# STOP EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+###########################################################################
+#
+# BELOW HERE, THINGS ARE DESIGNED NOT TO BE CHANGED
+#
+###########################################################################
+
+# We can either use interior regions, or breaklines
+
+if not interior_regions_data == []:
+    assert pt_areas is None, \
+        'Cannot define both ptAreas and non-empty interior regions'
+
+# bounding_polygon=su.read_polygon(bounding_polygon_file)
+
+bounding_polygon_and_tags = \
+    user_functions.read_boundary_tags_line_shapefile(
+        bounding_polygon_and_tags_file,
+        boundary_tags_attribute_name)
+bounding_polygon = bounding_polygon_and_tags[0]
+boundary_tags = bounding_polygon_and_tags[1]
+
+breaklines = su.readListOfBreakLines(breakline_files)
+(riverwalls, riverwall_par) = \
+    su.readListOfRiverWalls(riverwall_csv_files)  # {} # Dictionary
+
+if pt_areas is not None:
+    region_point_areas = su.readregion_point_areas(
+        pt_areas,
+        convert_length_to_area=True)
+else:
+    region_point_areas = None
+
+# Hack to override resolution
+# region_point_areas=[ region_point_areas[i][0:2]+[150*150*0.5] for i in \
+#                      range(len(region_point_areas))]
+
+# Redefine interior_regions to contain the polygon data + resolutions
+
+interior_regions = [[su.read_polygon(ir[0]), ir[1]] for ir in
+                    interior_regions_data]
+
+# Directory where sww files are written
+
+output_basedir = '../anuga/outputs/'
+
+# Directory where we store .txt versions of processed spatial inputs
+
+spatial_txt_output_dir = 'SpatialTxt/'
+
+# Deal with intersections in the bounding polygon / breaklines / riverwalls
+
+(bounding_polygon, breaklines, riverwalls) = \
+    su.add_intersections_to_domain_features(
+        bounding_polygon, breaklines, riverwalls,
+        point_movement_threshold=break_line_intersect_point_movement_threshold,
+        verbose=True)
+
+# Here we make a unique ID based on the all the mesh geometry inputs
+# This tells us if we need to regenerate partitions, or use old ones
+
+mesh_id_hash = hashlib.md5(json.dumps([
+    bounding_polygon,
+    interior_regions,
+    riverwalls,
+    breaklines,
+    region_point_areas,
+    default_res,
+    boundary_tags,
+])).hexdigest()
+
+# Fix the output tif bounding polygon
+
+if output_tif_bounding_polygon is None:
+    output_tif_bounding_polygon = bounding_polygon
+else:
+    output_tif_bounding_polygon = \
+        su.read_polygon(output_tif_bounding_polygon)
+
+# Make proj4string from projection information
+#
+
+if type(projection_information) is int:
+
+    # projection_information describes a UTM zone
+    # e.g. '+units=m +ellps=WGS84 +zone=47 +south=False +proj=utm '
+
+    if projection_information < 0:
+        proj4string = '+proj=utm +zone=' \
+            + str(abs(projection_information)) \
+            + ' +south +datum=WGS84 +units=m +no_defs'
+    else:
+        proj4string = '+proj=utm +zone=' + str(projection_information) \
+            + ' +datum=WGS84 +units=m +no_defs'
+elif type(projection_information) is str:
+    proj4string = projection_information
+else:
+    msg = 'Invalid projection information ' + \
+        ' --  must be a proj4string, or an integer' + \
+        ' defining a UTM zone [positive for northern hemisphere,' + \
+        ' negative for southern hemisphere]'
+    raise Exception(msg)
+
+# Set up directories etc
+
+if myid == 0:
+    runtime = time.strftime('%Y%m%d_%H%M%S', time.localtime())
+    for i in range(1, numprocs):
+        send(runtime, i)
+else:
+    runtime = receive(0)
+barrier()
+
+output_dir = output_basedir + 'RUN_' + str(runtime) + '_' + scenario
+partition_basedir = 'PARTITIONS/'
+partition_dir = partition_basedir + 'Mesh_' + str(mesh_id_hash)
+meshname = output_dir + '/mesh.tsh'
+
+
+def setup_directory_structure():
+    """
+    Make output directories, send important data to files
+    """
+
+    # if(myid==0):
+    # Make the base directory if it does not exist
+
+    try:
+        os.mkdir(output_basedir)
+    except:
+        pass
+
+    # Make the output directory
+
+    try:
+        os.mkdir(output_dir)
+    except:
+        pass
+
+    # Make the partitions directory (2 steps as it is recursive)
+
+    try:
+        os.mkdir(partition_basedir)
+    except:
+        pass
+
+    try:
+        os.mkdir(partition_dir)
+    except:
+        pass
+
+    # Make the spatialTxt directory
+
+    try:
+        os.mkdir(spatial_txt_output_dir)
+    except:
+        pass
+
+    # #########################################################################
+    # Export breaklines to txt for QC
+
+    def sav_lines(lines, filename='check.txt'):
+        new_a = numpy.array(lines)
+        numpy.savetxt(filename, new_a)
+        return
+
+    if myid == 0:
+        if interior_regions is not []:
+            for i in range(len(interior_regions)):
+                sav_lines(interior_regions[i][0],
+                          filename=spatial_txt_output_dir
+                          + 'interior_region_' + str(i) + '.txt')
+        if breaklines is not {}:
+            for i in range(len(breaklines)):
+                sav_lines(breaklines.values()[i],
+                          filename=spatial_txt_output_dir + 'breakline_'
+                          + str(i) + '.txt')
+        if riverwalls is not {}:
+            for i in range(len(riverwalls)):
+                sav_lines(riverwalls.values()[i],
+                          filename=spatial_txt_output_dir + 'riverwall_'
+                          + str(i) + '.txt')
+
+        sav_lines(bounding_polygon, filename=spatial_txt_output_dir
+                  + 'boundingPoly.txt')
+
+        if region_point_areas is not None:
+            sav_lines(region_point_areas,
+                      filename=spatial_txt_output_dir
+                      + 'regionPtAreas.txt')
+
+        # Copy all the txt representations of these files to the output
+        # directory
+
+        shutil.copytree(spatial_txt_output_dir, output_dir + '/'
+                        + spatial_txt_output_dir)
+    else:
+        pass
+
+    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/raster_outputs.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/raster_outputs.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/raster_outputs.py	(revision 9336)
@@ -0,0 +1,285 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+############################################################
+
+# Get the modules we need
+
+"""
+
+Code to make spatial outputs (GeoTifs)
+
+"""
+
+
+import sys
+import os
+import glob
+import numpy
+import gdal
+from pyproj import Proj
+from anuga.utilities import plot_utils as util
+from anuga.utilities import spatialInputUtil as su
+
+
+# @@@@@@@@@@@@@@@@@@#
+# STOP EDITING HERE
+# @@@@@@@@@@@@@@@@@@#
+#
+#
+
+def make_resampled_elevation(
+    elevation_raster,
+    raster_extent,
+    cell_size,
+    clip_polygon,
+    clip_polygon_layer,
+    proj4string,
+    tif_output_dir,
+):
+    """
+    Call gdalwarp to get the resampled elevation data
+    """
+
+    import platform
+    ostype = platform.system()
+
+    if ostype == 'Linux':
+        gdalwarp = 'gdalwarp'
+        gdal_rasterize = 'gdal_rasterize'
+    elif ostype == 'Windows':
+        gdalwarp = 'gdalwarp.exe'
+        gdalwarp = 'gdal_rasterize.exe'
+    else:
+        raise Exception('Unrecognized OS')
+
+    rast_xmin = raster_extent[0]
+    rast_xmax = raster_extent[1]
+    rast_ymin = raster_extent[2]
+    rast_ymax = raster_extent[3]
+    xres = cell_size
+    yres = cell_size
+
+    # allSRS='EPSG:'+str(EPSG_CODE)
+
+    all_srs = proj4string
+    new_dem = tif_output_dir + 'LIDAR_resampled.tif'
+    new_mask = tif_output_dir + 'Mask_resampled.tif'
+
+    extent_info = ' -te ' + str(rast_xmin) + ' ' + str(rast_ymin) + ' ' \
+        + str(rast_xmax) + ' ' + str(rast_ymax)
+    res_info = ' -tr ' + str(xres) + ' ' + str(yres)
+
+    # mkdir -p ElevationToPlot
+
+    gdalwarp_command = gdalwarp + ' -s_srs ' + all_srs + ' -t_srs ' \
+        + all_srs + extent_info + ' ' + res_info + ' ' \
+        + elevation_raster + ' ' + new_dem + ' -overwrite'
+
+    os.system(gdalwarp_command)
+
+    import shutil
+    shutil.copyfile(new_dem, new_mask)
+
+    gdal_rasterize_command = gdal_rasterize + ' -a_srs ' + all_srs \
+        + ' -burn 1 -l ' + clip_polygon_layer + ' ' + extent_info + ' ' \
+        + res_info + ' ' + clip_polygon + ' ' + new_mask
+    print gdal_rasterize_command
+    os.system(gdal_rasterize_command)
+
+    return
+
+
+#################################################################
+
+def gdal_calc_command(
+    stage,
+    depth,
+    elevation,
+    mask,
+    depth_threshold,
+):
+    """
+    Call gdal_calc to compute depth by draping
+    stage over the high res elevation data
+    """
+
+    gdal_calc = 'gdal_calc.py '
+    gdal_input = '-A ' + stage + ' -B ' + depth + ' -C ' + elevation \
+        + ' -D ' + mask
+    calc_command = " --calc '(A-C)*(A>C)*(B>" + str(depth_threshold) \
+        + ")*D' "
+
+    # rastOut=$swwTifDir$swwBase'_depth_high_res.tif'
+
+    rast_out = os.path.dirname(depth) + '/HighRes_' \
+        + os.path.basename(depth)
+
+    gdalcalc_command = gdal_calc + gdal_input + calc_command \
+        + '--outfile ' + rast_out
+    print gdalcalc_command
+    os.system(gdalcalc_command)
+
+    return
+
+
+############################################################
+
+def make_me_some_tifs(
+    sww_file,
+    bounding_polygon,
+    proj4string,
+    my_time_step='collected_max',
+    tif_output_subdir='/TIFS/',
+    cell_size=5.0,
+    make_highres_drape_plot=False,
+    elevation_raster=None,
+    depth_threshold=None,
+    clip_polygon=None,
+    clip_polygon_layer=None,
+):
+    """
+
+    ### INPUT DATA
+    - swwFile -- Full path name of sww to read outputs from
+    - bounding_polygon -- ANUGA's bounding polygon, or another polygon to clip
+      the rasters to
+    - proj4string defining the coordinate system
+    - myTimeStep -- option from util.Make_Geotif
+      use 'max' to plot the maxima
+      use [0, 5, 10] to plot the
+      use 'collected_max' to read the csv outputs from
+        collect_max_quantities_operator
+    - tif_outputdir -- Write outputs to this folder inside the swwFile's
+      directory (MUST INCLUDE TRAILING SLASH /)
+    - CellSize -- Desired raster cellSize (m)
+    - make_highres_drape_plot -- True/False, Make a high-res drape plot?
+    - elevation_raster -- Filename of elevation raster for 'high-res-drape'
+      depth plot [from subtracting stage from high res topography]
+    - depth_threshold -- Depth threshold for high-res-drape' plot (m). If
+      ANUGA's triangle has depth
+      < depth_threshold, then depth=0 in all high-res cells in that triangle
+    - clipPolygon -- Polygon to clip 'high-res-drape' plot. Must be
+      provided if make_highres_drape_plot==True (can use bounding polygon or
+      another choice)
+    - clipPolygonLayer -- Layer for above as used by gdal (usually shapefile
+      name without .shp)
+
+    ## OUTPUT
+    Nothing is returned, but tifs are made in tif_output_subdir inside the
+    swwFile directory
+
+    """
+
+    # Convert utm_zone to proj4string
+    # p=Proj(proj='utm', south=(utm_zone<0.),
+    #        zone=abs(utm_zone), ellps='WGS84')
+    # proj4string = p.srs
+
+    tif_output_dir = os.path.dirname(sww_file) + tif_output_subdir
+
+    # Make the geotifs
+
+    if my_time_step == 'collected_max':
+
+        # Read max quantity output files inputs
+
+        max_qfiles = glob.glob(os.path.dirname(sww_file)
+                               + '/*_UH_MAX.csv')
+        if len(max_qfiles) == 0:
+            raise Exception(
+                'Cannot find any files containing collected maxima')
+        for i in range(len(max_qfiles)):
+            if i == 0:
+                max_qs = numpy.genfromtxt(max_qfiles[i], delimiter=',')
+            else:
+                extra_data = numpy.genfromtxt(max_qfiles[i], delimiter=',')
+                max_qs = \
+                    numpy.vstack([max_qs, extra_data])
+
+        # Make the geotiff's
+
+        for (i, quant) in enumerate(['stage', 'depth', 'velocity',
+                                     'depthIntegratedVelocity']):
+
+            # FIXME: The interpolation function is remade for every quantity,
+            # since only a 3 column array can be passed to Make_Geotif Could
+            # make it fast (do it only once) by changing Make_Geotif
+
+            tmp_arr = max_qs[:, [0, 1, i + 2]]
+            util.Make_Geotif(
+                tmp_arr,
+                output_quantities=[quant + '_MAX'],
+                cell_size=cell_size,
+                proj4string=proj4string,
+                verbose=True,
+                bounding_polygon=bounding_polygon,
+                output_dir=tif_output_dir,
+            )
+
+        # Also plot elevation + friction
+
+        util.Make_Geotif(
+            sww_file,
+            my_time_step=[0],
+            output_quantities=['elevation', 'friction'],
+            cell_size=cell_size,
+            proj4string=proj4string,
+            verbose=True,
+            bounding_polygon=bounding_polygon,
+            output_dir=tif_output_dir,
+        )
+    else:
+        util.Make_Geotif(
+            sww_file,
+            my_time_step=my_time_step,
+            output_quantities=[
+                'depth',
+                'stage',
+                'elevation',
+                'velocity',
+                'depthIntegratedVelocity',
+                'friction',
+            ],
+            cell_size=cell_size,
+            proj4string=proj4string,
+            verbose=True,
+            bounding_polygon=bounding_polygon,
+            output_dir=tif_output_dir,
+        )
+
+    # Early finish
+
+    if not make_highres_drape_plot:
+        return
+
+    # Get extent of geotifs
+
+    sample_rast = glob.glob(tif_output_dir + '*.tif')[0]
+    raster_extent = su.getRasterExtent(sample_rast)
+
+    # Make the resampled elevation data
+
+    make_resampled_elevation(
+        elevation_raster,
+        raster_extent,
+        cell_size,
+        clip_polygon,
+        clip_polygon_layer,
+        proj4string,
+        tif_output_dir,
+    )
+
+    elevation = glob.glob(tif_output_dir + 'LIDAR_resampled*')[0]
+    mask = glob.glob(tif_output_dir + 'Mask_resampled*')[0]
+    if my_time_step == 'collected_max':
+        depth = glob.glob(tif_output_dir + 'PointData_depth_*')[0]
+        stage = glob.glob(tif_output_dir + 'PointData_stage_*')[0]
+    else:
+        depth = glob.glob(tif_output_dir + '*_depth_max.tif')[0]
+        stage = glob.glob(tif_output_dir + '*_stage_max.tif')[0]
+
+    # Call gdal_calc
+
+    gdal_calc_command(stage, depth, elevation, mask, depth_threshold)
+
+    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/run_model.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/run_model.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/run_model.py	(revision 9336)
@@ -0,0 +1,180 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+# ------------------------------------------------------------------------------
+# Import necessary modules
+# ------------------------------------------------------------------------------
+# Standard modules
+
+"""
+
+Script for running a semi-generic anuga model
+
+It's supposed to be organised in a fairly conceptual manner -- a collection of
+clearly named single-line function calls which do the main part of the work. To
+customize it, most edits will be required in the associated setup_* scripts
+
+If you want to add/adjust things , consider putting detailed code as functions
+in user_functions.py , and calling them from here -- to keep it clean
+
+Gareth Davies
+Geoscience Australia, 2014-present
+
+"""
+
+
+import os
+import time
+
+# ANUGA modules
+
+import anuga
+from anuga_parallel import myid, numprocs, finalize, barrier
+from anuga.operators.collect_max_quantities_operator import \
+    collect_max_quantities_operator
+
+# Routines to setup the domain
+
+import setup_boundary_conditions
+import setup_rainfall
+import setup_inlets
+import setup_mesh
+import setup_initial_conditions
+
+# Routines defined by the user
+
+import user_functions
+
+# Functions to make geotifs
+
+import raster_outputs
+
+# Record the time so we can report how long the simulation takes
+
+t0 = time.time()
+
+# Get key data for the simulation
+
+import project
+
+###############################################################################
+#
+# SETUP DOMAIN
+#
+###############################################################################
+
+print 'Making domain'
+domain = setup_mesh.setup_mesh(project, setup_boundary_conditions)
+
+##########################################################################
+#
+# SETUP INITIAL CONDITIONS
+#
+##########################################################################
+
+print 'Making initial conditions'
+setup_initial_conditions.setup_initial_conditions(domain, project)
+
+##########################################################################
+#
+# SETUP FORCING TERMS
+#
+##########################################################################
+# print 'Making rainfall conditions'
+# setup_rainfall.setup_rainfall(domain, project)
+# print 'Making inlets '
+# setup_inlets.setup_inlets(domain, project)
+
+##########################################################################
+#
+# SETUP BOUNDARY CONDITIONS
+#
+##########################################################################
+
+print 'Making boundary conditions '
+setup_boundary_conditions.setup_boundary_conditions(domain, project)
+
+##########################################################################
+#
+# STORE MAX QUANTITIES
+#
+##########################################################################
+
+max_quantities = collect_max_quantities_operator(
+    domain,
+    update_frequency=project.max_quantity_update_frequency,
+    collection_start_time=project.max_quantity_collection_start_time,
+    # Set stored velocities to zero if height<velocity_zero_height
+    velocity_zero_height=1.0e-03)
+
+##########################################################################
+#
+# EVOLVE IN TIME
+#
+##########################################################################
+
+print 'Evolving'
+
+barrier()
+for t in domain.evolve(yieldstep=project.yieldstep,
+                       finaltime=project.finaltime):
+    if myid == 0:
+        print domain.timestepping_statistics()
+
+    # ##########################################################################
+    # Print various things (optional)
+    #
+    # Compute the volume
+
+    domain.report_water_volume_statistics()
+
+    # Report cells with a small local timestep
+    # domain.report_cells_with_small_local_timestep()
+
+    # Local function to print the velocity statistics
+    # user_functions.print_velocity_statistics(domain,MaxQuantities)
+
+barrier()
+
+###############################################################################
+#
+# POST-SIMULATION PROCESSING
+#
+###############################################################################
+
+# Write max quantities to a file
+
+max_quantity_file_start = domain.get_datadir() + '/Max_quantities_'
+max_quantities.export_max_quantities_to_csv(max_quantity_file_start)
+
+# Merge the parallel SWW files together
+
+os.chdir(project.output_dir)
+if myid == 0 and numprocs > 1:
+    print 'Number of processors %g ' % numprocs
+    print 'That took %.2f seconds' % (time.time() - t0)
+    print 'Communication time %.2f seconds' % domain.communication_time
+    print 'Reduction Communication time %.2f seconds' \
+        % domain.communication_reduce_time
+    print 'Broadcast time %.2f seconds' \
+        % domain.communication_broadcast_time
+
+    anuga.utilities.sww_merge.sww_merge_parallel(
+        project.scenario,
+        np=numprocs, verbose=True, delete_old=True)
+
+# FIXME: Make geotif rasters
+
+if myid == 0:
+    try:
+        raster_outputs.make_me_some_tifs(
+            sww_file='./' + project.scenario + '.sww',
+            bounding_polygon=project.bounding_polygon,
+            proj4string=project.proj4string,
+            cell_size=project.output_tif_cellsize)
+    except:
+        print 'GeoTif creation failed, you can try manually using' + \
+              ' raster_outputs.py or anuga.utilities.plot_utils.Make_Geotif'
+
+barrier()
+finalize()
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_boundary_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_boundary_conditions.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_boundary_conditions.py	(revision 9336)
@@ -0,0 +1,74 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+
+THE USER MUST EDIT THIS FILE TO SET THE BOUNDARY CONDITIONS
+
+"""
+
+
+import numpy
+import scipy
+import scipy.interpolate
+import anuga
+import anuga.utilities.spatialInputUtil as su
+
+
+def setup_boundary_conditions(domain, project):
+    """
+    THE USER MUST MODIFY THIS TO SETUP THE BOUNDARIES AS DESIRED
+
+    1) First make sure the required boundary conditions are created
+    (see examples below, modify or add new ones as required)
+
+    2) Next make sure the 'boundary_tags_and_conditions' dictionary
+    associates the correct 'tags' (contained as attribute names in the
+    bounding_polygon_and_tags_file) with the correct boundary condition type
+    """
+
+# @@@@@@@@@@@@@@@@@#
+# START EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+    # Example:
+    # "Standard" boundary conditions
+
+    br = anuga.Reflective_boundary(domain)
+
+    # Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Be
+    # careful with this one
+
+    # Example:
+    # Transmissive with stage being set by a function 'stage_function'
+
+    def stage_function(t):
+        return min(max(t - 60., 0.) / 60., 1.)
+
+    import anuga.shallow_water.boundaries as asb
+    b_t_stage_function = \
+        asb.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(
+            domain,
+            stage_function)
+
+    # Here we associate each tag in boundary_info with its boundary condition
+    # defined above
+
+    boundary_tags_and_conditions = {'Ocean': b_t_stage_function,
+                                    'Br': br}
+
+# @@@@@@@@@@@@@@@@@#
+# STOP EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+    # Check that all tags in boundary_info have been set
+
+    for i in range(len(project.boundary_tags)):
+        tag = project.boundary_tags.keys()[i]
+        if not (tag in boundary_tags_and_conditions.keys()):
+            msg = 'Need to set boundary_tags_and_conditions for tag = ' \
+                + tag
+            raise Exception(msg)
+
+    domain.set_boundary(boundary_tags_and_conditions)
+
+    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_initial_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_initial_conditions.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_initial_conditions.py	(revision 9336)
@@ -0,0 +1,276 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+"""
+
+THIS SHOULD BE CHANGED BY THE USER TO SET THE ELEVATION, FRICTION, INITIAL
+STAGE etc
+
+These settings are defined in the function setup_initial_conditions
+
+"""
+
+
+import numpy
+import glob
+import anuga
+from anuga.utilities import quantity_setting_functions as qs
+from anuga.utilities import spatialInputUtil as su
+
+
+def setup_initial_conditions(domain, project):
+    """
+    Set the initial conditions for ANUGA
+
+    INPUT: domain  = the anuga domain
+    project = the project module
+
+    """
+
+    # #########################################################################
+    # INSTRUCTIONS ON USAGE
+    #
+    # Here we set the initial quantities, optionally applying different types
+    # of interpolation in different polygons.
+    #
+    # To do this, each quantity needs an ordered list of [polygon, value] pairs
+    #  (where we will set xy points in 'polygon' using 'value').
+    #
+    # Possible values of 'polygon' are:
+    #  1) An ANUGA polygon
+    #  2) A filename with .shp extension, or a csv file in anuga_polygon format
+    #  (with any extension)
+    #  3) None -- apply nowhere
+    #  4) 'Extent' (only valid if value is a raster filename) -- apply over the
+    #  full raster extent
+    #  5) 'All' -- apply everywhere -- must be the last polygon-value pair
+    #
+    # Possible values of 'value' are:
+    # 1) A constant
+    # 2) A 3-column numpy array with x,y,z (then nearest neighbour
+    # interpolation is used)
+    # 3) A gdal-compatible raster (pixel-lookup is used)
+    # 4) Any function f(x,y). Note the function must act on 'ANUGA-transformed'
+    # xy coordinates which have
+    #    min(x), min(y) = (0.,0.)
+    #
+    # If polygons overlap, ones earlier in the list take preference
+    # e.g.
+    #   [ [polygon_1, value_1], [polygon_2, value_2],
+    #     [polygon_3, value_3], ['All', value_4] ]
+    # will set all points in polygon_1 using value_1, then all points in
+    # polygon_2 which were not in polygon_1 using value_2, then all points in
+    # polygon_3 which were not previously set using value_3, then everything
+    # remaining with value_4
+    #
+    #
+    # For more information, see (in anuga.utilities.quantity_setting_functions)
+    #   composite_quantity_setting_function
+    #
+
+    # @@@@@@@@@@@@@@@@@#
+    # START EDITING HERE
+    # @@@@@@@@@@@@@@@@@#
+
+    # #########################################################################
+    #
+    # PRELIMINARY DEFINITIONS
+    #
+    # #########################################################################
+
+    # If the earthquake slip is given then convert to deformation here
+    # So that we can add it to both elevation AND stage
+
+    make_earthquake_deformation_from_source_model = False
+
+    if make_earthquake_deformation_from_source_model:
+        source_model = '../anuga/source/XXXXX.csv'
+
+        # Make an earthquake source model
+
+        eq_function = okada_tsunami.earthquake_source(
+            source=None,
+            filename=source_model,
+            domain=domain,
+            lon_lat_degrees=True,
+            lon_before_lat=False,
+            proj4string=project.proj4string,
+            verbose=False,
+        )
+    else:
+
+        # Make a function that is zero everywhere
+        # If we pass this to add_quantity, nothing should happen
+
+        def eq_function(x, y):
+            return x * 0.
+
+    # #########################################################################
+    #
+    # Set elevation.
+    #
+    # #########################################################################
+    # The 'Patong 10m' dataset
+
+    elevation_patong_10m = \
+        numpy.genfromtxt(
+            '../anuga/topographies/patong_10m_grid_msl_Project.txt',
+            delimiter=',', skip_header=1)
+
+    # The 'saddle 10m' dataset
+
+    sf = '../anuga/topographies/patong_10m_small_grid_for_anuga_sub_Project.txt'
+    elevation_saddle_10m = \
+        numpy.genfromtxt(sf, delimiter=',', skip_header=1)
+
+    # The '1s grid' dataset
+
+    elevation_patong_1s = \
+        numpy.genfromtxt(
+            '../anuga/topographies/patong_bay_1s_grid_Project.txt',
+            delimiter=',', skip_header=1)
+
+    # The background data
+
+    elevation_default = \
+        numpy.genfromtxt(
+            '../anuga/topographies/andaman_3s_grid_Clip2_Project.txt',
+            delimiter=',', skip_header=1)
+
+    # Set_quantity data
+
+    elevation_data = [['../anuga/polygons/patong_10m.txt',
+                       elevation_patong_10m],
+                      ['../anuga/polygons/saddle_10m.txt',
+                       elevation_saddle_10m],
+                      ['../anuga/polygons/patong_1s.txt',
+                       elevation_patong_1s], ['All', elevation_default]]
+
+    # Define a function for add_quantity in the same way
+    # [ ['All', 0.] ] will do nothing
+
+    elevation_additions = [['All', eq_function]]
+
+    # #############################################################################
+    #
+    # Set friction
+    #
+    # #############################################################################
+
+    default_friction = 0.025
+
+    # Set_quantity data
+
+    friction_data = [['All', default_friction]]
+
+    # Add_quantity data
+
+    friction_additions = [['All', 0.]]
+
+    # ##############################################################################
+    #
+    # Set stage
+    #
+    # ##############################################################################
+
+    # Set_quantity data
+
+    stage_data = [['All', project.tide]]
+
+    # Add_quantity data
+
+    stage_additions = [['All', eq_function]]
+
+    # ##############################################################################
+    #
+    # Set depth integrated velocity (called "momentum" in ANUGA although that
+    # is incorrect")
+    #
+    # ##############################################################################
+
+    # Set_quantity data
+
+    xmomentum_data = [['All', 0.]]
+
+    # Add_quantity data
+
+    xmomentum_additions = [['All', 0.]]
+
+    # Set_quantity data
+
+    ymomentum_data = [['All', 0.]]
+
+    # Add_quantity data
+
+    ymomentum_additions = [['All', 0.]]
+
+    # @@@@@@@@@@@@@@@@@#
+    # STOP EDITING HERE
+    # @@@@@@@@@@@@@@@@@#
+
+    # ###############################################################################
+    #
+    # Set the quantities here
+    #
+    # ###############################################################################
+
+    # To allow for 'discontinuous' elevation or other quantities, set
+    # location='centroids'
+
+    elevation_function = \
+        qs.composite_quantity_setting_function(elevation_data, domain)
+    domain.set_quantity('elevation', elevation_function,
+                        location='vertices')
+    elevation_addition_function = \
+        qs.composite_quantity_setting_function(elevation_additions,
+                                               domain)
+    domain.add_quantity('elevation', elevation_addition_function)
+
+    friction_function = \
+        qs.composite_quantity_setting_function(friction_data, domain)
+    domain.set_quantity('friction', friction_function,
+                        location='centroids')
+    friction_addition_function = \
+        qs.composite_quantity_setting_function(friction_additions,
+                                               domain)
+    domain.add_quantity('friction', friction_addition_function)
+
+    stage_function = qs.composite_quantity_setting_function(stage_data,
+                                                            domain)
+    domain.set_quantity('stage', stage_function, location='centroids')
+    stage_addition_function = \
+        qs.composite_quantity_setting_function(stage_additions, domain)
+    domain.add_quantity('stage', stage_addition_function)
+
+    xmomentum_function = \
+        qs.composite_quantity_setting_function(xmomentum_data, domain)
+    domain.set_quantity('xmomentum', xmomentum_function,
+                        location='centroids')
+    xmomentum_addition_function = \
+        qs.composite_quantity_setting_function(xmomentum_additions,
+                                               domain)
+    domain.add_quantity('xmomentum', xmomentum_addition_function)
+
+    ymomentum_function = \
+        qs.composite_quantity_setting_function(ymomentum_data, domain)
+    domain.set_quantity('ymomentum', ymomentum_function,
+                        location='centroids')
+    ymomentum_addition_function = \
+        qs.composite_quantity_setting_function(ymomentum_additions,
+                                               domain)
+    domain.add_quantity('ymomentum', ymomentum_addition_function)
+
+    # #########################################################################
+    #
+    # Add Riverwalls [ must happen after distribute(domain) in parallel ]
+    #
+    # #########################################################################
+
+    if not project.riverwalls == {}:
+        domain.riverwallData.create_riverwalls(project.riverwalls,
+                                               project.riverwall_par)
+        domain.riverwallData.export_riverwalls_to_text(
+            output_dir=project.output_dir + '/' +
+            project.spatial_txt_output_dir)
+
+    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_inlets.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_inlets.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_inlets.py	(revision 9336)
@@ -0,0 +1,69 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+
+THE USER CAN ADJUST THIS TO APPLY INLET INFLOWS
+
+Set up inflows
+
+Read a discharge input timeseries, and impose it
+
+"""
+
+
+import numpy
+import anuga
+from anuga.utilities import spatialInputUtil as su
+import scipy
+import scipy.interpolate
+
+# @@@@@@@@@@@@@@@@@#
+# START EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+##########################################################################
+##
+# LOWER MANGAHAN INFLOWS
+# LowerMangahanDischarge=numpy.genfromtxt('../../Boundary_timeseries/Lower_Marikina_Ondoy.csv',
+#                                        delimiter=",",skip_header=1)
+# Only keep time and discharge. Time starts from noon, 2005-09-25
+# LowerMangahanDischarge=LowerMangahanDischarge[:,[0,2]]
+# Convert to function
+# LowerMangahan_QFun=scipy.interpolate.interp1d(LowerMangahanDischarge[:,0],
+#                                        LowerMangahanDischarge[:,1],kind='linear')
+# LowerMangahan_line=su.readShp_1LineGeo('../../GIS/Inlets/LowerMangahan.shp')
+
+###############################################################################
+#
+# This determines the discharge inflows
+
+# inlet_discharges=[
+#                  [LowerMangahan_line, LowerMangahan_QFun]
+#                   ]
+
+inlet_discharges = []
+
+
+# @@@@@@@@@@@@@@@@@#
+# STOP EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+###############################################################################
+#
+# Below here nothing should change
+#
+###############################################################################
+
+def setup_inlets(domain, project):
+    """
+    Add inlets to domain
+    """
+
+    for i in range(len(inlet_discharges)):
+        line = inlet_discharges[i][0]
+        qfun = inlet_discharges[i][1]
+
+        # Add inlet
+
+        anuga.Inlet_operator(domain, line, qfun)
+    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_mesh.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_mesh.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_mesh.py	(revision 9336)
@@ -0,0 +1,256 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+"""
+
+THIS IS NOT DESIGNED TO BE CHANGED BY THE USER
+
+Controlling variables are defined in
+
+project.py (polygons, resolutions etc)
+setup_boundary_conditions.py (Data used to make boundary tags)
+
+
+"""
+
+
+import glob
+import os
+from os.path import join
+import gc
+
+import anuga
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+from anuga.utilities import spatialInputUtil as su
+from anuga_parallel.parallel_api import pypar_available
+if pypar_available:
+    from anuga import sequential_distribute_load
+    from anuga import sequential_distribute_dump
+    from anuga_parallel.sequential_distribute import \
+        sequential_distribute_load_pickle_file
+
+# @@@@@@@@@@@@@@@@@#
+# STOP EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+#
+# Yes, I mean stop - in theory it's not required.
+#
+
+verbose = True
+
+
+def build_mesh(project, setup_boundary_conditions):
+    """
+
+    This is executed by processor 0 to build the mesh.
+
+    """
+
+    # Create the boundary tags
+    # boundary_tags=\
+    #    make_boundary_tags(project.bounding_polygon,
+    #                       setup_boundary_conditions)
+
+    # Ensure mesh_breaklines include riverwalls and breaklines
+
+    mesh_breaklines = \
+        su.combine_breakLines_and_riverWalls_for_mesh(project.breaklines,
+                                                      project.riverwalls)
+
+    # Make the mesh
+
+    anuga.create_mesh_from_regions(
+        project.bounding_polygon,
+        boundary_tags=project.boundary_tags,
+        maximum_triangle_area=project.default_res,
+        filename=project.meshname,
+        interior_regions=project.interior_regions,
+        use_cache=False,
+        verbose=verbose,
+        breaklines=mesh_breaklines,
+        regionPtArea=project.region_point_areas,
+    )
+
+    # Make the domain using the mesh
+
+    domain = anuga.create_domain_from_file(project.meshname)
+
+    # Print some stats about mesh and domain
+
+    print 'Number of triangles = ', len(domain)
+    print 'The extent is ', domain.get_extent()
+    print domain.statistics()
+
+    # Print info on the smallest triangles
+
+    small_areas = domain.areas.argsort()
+    print ''
+    print 'LOCATIONS OF TRIANGLES WITH SMALLEST AREAS'
+    for i in range(10):
+        j = small_areas[i]
+        x = domain.centroid_coordinates[j, 0] \
+            + domain.geo_reference.xllcorner
+        y = domain.centroid_coordinates[j, 1] \
+            + domain.geo_reference.yllcorner
+        print '  Area ' + str(domain.areas[j]) + ' location: ' \
+            + str(round(x, 1)) + ',' + str(round(y, 1))
+    print ''
+
+    return domain
+
+
+##########################################################################
+
+def setup_mesh(project, setup_boundary_conditions):
+    """
+    Code to make the mesh (initial domain)
+
+    The geometry is made on processor 0, then dumped and reloaded
+    This reduces the memory demands
+
+    INPUT: project == the project module
+
+    OUTPUT: domain
+    """
+
+    if myid == 0:
+
+        # Make sure files/folders to store everything are there
+
+        project.setup_directory_structure()
+
+        # Copy code files to output dir. This should be 'nearly the first'
+        # thing we do
+
+        model_files = glob.glob('*.py') + glob.glob('*.sh')
+        for model_file in model_files:
+            anuga.copy_code_files(project.output_dir, model_file)
+
+        if verbose:
+            print 'Hello from processor ', myid
+
+        #
+        # HERE, WE MAKE/PARTITION/READ THE MESH
+        # This can lead to memory savings in parallel runs
+        # (possibly because of limitations of python memory management)
+        #
+
+        # Let's see if we have already pickled this domain
+
+        pickle_name = 'domain' + '_P%g_%g.pickle' % (1, 0)
+        pickle_name = join(project.partition_dir, pickle_name)
+
+        if os.path.exists(pickle_name):
+            if verbose:
+                print 'Saved domain seems to already exist'
+        else:
+            if verbose:
+                print 'CREATING PARTITIONED DOMAIN'
+            domain = build_mesh(project, setup_boundary_conditions)
+
+            # If pypar is not available don't use sequential_distribute stuff
+            # (it will fail)
+
+            if pypar_available:
+                if verbose:
+                    print 'Saving Domain'
+                sequential_distribute_dump(domain, 1,
+                                           partition_dir=project.partition_dir,
+                                           verbose=verbose)
+
+        # If pypar is not available don't use sequential_distribute stuff (it
+        # will fail)
+
+        if pypar_available:
+
+            # Now partition the domain
+
+            par_pickle_name = 'domain' + '_P%g_%g.pickle' % (numprocs,
+                                                             0)
+            par_pickle_name = join(project.partition_dir,
+                                   par_pickle_name)
+            if os.path.exists(par_pickle_name):
+                if verbose:
+                    print 'Saved partitioned domain seems to already exist'
+            else:
+                if verbose:
+                    print 'Load in saved sequential pickled domain'
+                domain = \
+                    sequential_distribute_load_pickle_file(
+                        pickle_name, np=1, verbose=verbose)
+
+                if verbose:
+                    print 'Dump partitioned domains'
+                sequential_distribute_dump(
+                    domain, numprocs,
+                    partition_dir=project.partition_dir, verbose=verbose)
+
+            # This can reduce the memory demands if the domain was made above
+            # Even better to have an existing partition (i.e. stop here and
+            # rerun)
+
+            domain = None
+            gc.collect()
+    else:
+
+        domain = None
+        if verbose:
+            print 'Hello from processor ', myid
+
+    barrier()
+
+    # print 'Distributing domain'
+    # domain=distribute(domain)
+    # barrier()
+
+    # If pypar is not available don't use sequential_distribute stuff (it will
+    # fail)
+
+    if pypar_available:
+        if myid == 0:
+            print 'LOADING PARTITIONED DOMAIN'
+
+        domain = \
+            sequential_distribute_load(
+                filename=join(project.partition_dir, 'domain'),
+                verbose=verbose)
+
+    # #########################################################################
+    # Set output directories
+    # #########################################################################
+
+    domain.set_name(project.scenario)  # Name of sww file
+    domain.set_datadir(project.output_dir)  # Store sww output here
+
+    # Needs more changes for this to work
+    # domain.set_checkpointing(checkpoint_time=project.checkpoint_time)
+
+    # #########################################################################
+    # Miscellanious numerics
+    # #########################################################################
+
+    domain.set_flow_algorithm(project.flow_algorithm)
+
+    # Force zero beta values [hopefully integrated into source]
+    # print 'Warning: Forcing everything to first order'
+    # domain.beta_w=0.
+    # domain.beta_uh=0.
+    # domain.beta_vh=0.
+    # domain.beta_w_dry=0.
+    # domain.beta_uh_dry=0.
+    # domain.beta_vh_dry=0.
+
+    # Adjust velocity computation for max quantities
+    # domain.velocity_protection=1.0e-05
+
+    # Adjust CFL
+    # domain.set_CFL(0.9)
+
+    # Optionally store vertex values uniquely (large file sizes!)
+
+    domain.set_store_vertices_uniquely(project.store_vertices_uniquely)
+
+    if project.use_local_extrapolation_and_flux_updating:
+        domain.set_local_extrapolation_and_flux_updating()
+
+    return domain
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_rainfall.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_rainfall.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/setup_rainfall.py	(revision 9336)
@@ -0,0 +1,80 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+
+THE USER CAN ADJUST THIS TO APPLY RAINFALL
+
+Set rainfall
+
+Make two-column array with Time (s), rain (m/s) , and pass to linear
+interpolation function
+
+"""
+
+
+import numpy
+import anuga
+
+# @@@@@@@@@@@@@@@@@#
+# START EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+# idealised_constant_rain=False
+# if(idealised_constant_rain):
+# Good for testing / finding slowpoints in mesh / etc
+#    #
+# This series has:
+# 90mm of rain per hour initially, then nothing, 24 hours in total
+# (note 24 hours requires 25 interpolation points!)
+# rain_rate=90. #mm/hr
+# rain_rate=5. #mm/hr
+#    rain_timeseries=numpy.array(
+#        [numpy.linspace(0., 24*3600.,25),\
+#         10*[rain_rate/(3600.*1000.)]+15*[0.]]).transpose()
+# else:
+# Use data
+#    #
+#    rain_timeseries=numpy.genfromtxt(
+#        '../../Boundary_timeseries/science_garden_ondoy_rain.csv',
+#                                     skip_header=2,
+#                                     delimiter=',')
+# Add time, beginning 8.30 am 2009-09-25
+#    rain_timeseries=numpy.vstack(
+#        [numpy.linspace(0., (rain_timeseries.shape[0]-1)*3600.,
+#                             rain_timeseries.shape[0]),
+#                             rain_timeseries[:,5]]).transpose()
+# Adjust time and rainfall units
+# rain_timeseries[:,0]=rain_timeseries[:,0]-3.5*3600. # zero @ noon 25-09-2009
+# rain_timeseries[:,1]=rain_timeseries[:,1]/(3600.*1000) # mm/hr --> m/s
+#
+# print 'Reducing rainfall by factor of 5'
+# rain_timeseries[:,1]=rain_timeseries[:,1]/5.
+#
+#
+# Create function from timeseries with nearest-neighbour interpolation
+# Note that linear interpolation will not satisfy the
+# 'total-rainfall-over-the-hour' constraint
+# implied by the data
+# import scipy.interpolate
+# myrain=scipy.interpolate.interp1d(
+#    rain_timeseries[:,0], rain_timeseries[:,1],\
+#    kind='nearest')
+
+myrain = 0.
+
+
+# @@@@@@@@@@@@@@@@@#
+# STOP EDITING HERE
+# @@@@@@@@@@@@@@@@@#
+
+def setup_rainfall(domain, project):
+    """
+    Function to add rainfall operators to the domain
+    """
+
+    if not myrain == 0.:
+        rainin_term = \
+            anuga.operators.rate_operators.Rate_operator(domain,
+                                                         rate=myrain)
+
+    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/user_functions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/user_functions.py	(revision 9336)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami2/scripts/user_functions.py	(revision 9336)
@@ -0,0 +1,134 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+"""
+
+Home for various 'user-defined' functions.
+
+The user can define any new functions they need in this script, and make
+small modifications to the main scripts to use it.
+
+This is usually nicer than making detailed hacky changes to the main scripts
+
+If you have a nice function, consider pushing it to the ANUGA source
+
+"""
+
+
+import anuga
+from anuga_parallel import myid, numprocs, finalize, barrier
+
+
+def print_velocity_statistics(domain, max_quantities):
+    """
+    Utility for printing velocity stats in the evolve loop
+    """
+
+    # Print velocity statistics
+
+    for i in range(numprocs):
+        if myid == i:
+
+            # Peak speed info
+
+            xx = domain.quantities['xmomentum'].centroid_values
+            yy = domain.quantities['ymomentum'].centroid_values
+            dd = domain.quantities['stage'].centroid_values \
+                - domain.quantities['elevation'].centroid_values
+            dd = dd * (dd > 1.0e-03) + 1.0e-03 * (dd <= 1.0e-03)
+            vv = 1 / dd * (xx ** 2 + yy ** 2) ** 0.5
+            vv = vv * (dd > 1.0e-03)
+            print '    Processor ', myid
+            print '    @ Peak velocity is: ', vv.max(), vv.argmax()
+            print '     &- MaxSpeedHistory: ', \
+                max_quantities.max_speed.max()
+            print '     %- FUF: ', domain.flux_update_frequency.mean()
+        else:
+            pass
+        barrier()
+
+    # Make a newline
+
+    if myid == 0:
+        print ''
+
+
+############################################################
+
+def read_boundary_tags_line_shapefile(shapefile_name,
+                                      tag_attribute='bndryTag'):
+    """
+    Read in the boundary lines + tags from an appropriately structured line
+    shapefile
+
+    Return the bounding polygon + dict with tags + indices
+    """
+
+    # Step 1: Read the data
+
+    import ogr
+
+    data_source = ogr.Open(shapefile_name)
+
+    assert data_source.GetLayerCount() == 1, \
+        ' Bounding polygon + boundary tags file can only have 1 layer'
+
+    layer = data_source.GetLayer()
+
+    boundary_line_and_tags = []
+
+    for feature in layer:
+        line = feature.GetGeometryRef().GetPoints()
+        line = [list(l) for l in line]
+        tag = feature.GetField(tag_attribute)
+        boundary_line_and_tags.append([line, tag])
+
+    # Step 2: Convert to bounding polygon + boundary tags in ANUGA format
+
+    boundary_segnum = len(boundary_line_and_tags)
+
+    # Initial values
+
+    bounding_polygon = boundary_line_and_tags[0][0]
+    boundary_tags = \
+        {boundary_line_and_tags[0][1]: range(len(bounding_polygon) - 1)}
+
+    #
+    # import pdb
+    # pdb.set_trace()
+
+    for i in range(boundary_segnum - 1):
+
+        # Find the line which starts at the end point of bounding_polygon
+
+        for j in range(boundary_segnum):
+            blj = boundary_line_and_tags[j][0]
+            if blj[0] == bounding_polygon[-1]:
+
+                # Append all but the first point of blj to bounding_polygon
+
+                old_lbp = len(bounding_polygon)
+                extra_polygon = [blj[k] for k in range(1, len(blj))]
+                bounding_polygon = bounding_polygon + extra_polygon
+
+                # If we are on the last segment, then drop the final point
+                # (for consistency with ANUGA format)
+
+                if i == boundary_segnum - 2:
+                    bounding_polygon = [bounding_polygon[k] for k in
+                                        range(len(bounding_polygon) - 1)]
+                    new_tags = range(old_lbp - 1, len(bounding_polygon))
+                else:
+                    new_tags = range(old_lbp - 1, len(bounding_polygon)
+                                     - 1)
+
+                # Either the key already exists and we append edge indices, OR
+                # we add the key and the edge indices
+
+                key = boundary_line_and_tags[j][1]
+                if (key in boundary_tags.keys()):
+                    boundary_tags[key] = boundary_tags[key] + new_tags
+                else:
+                    boundary_tags[key] = new_tags
+
+    return (bounding_polygon, boundary_tags)
