Index: uga_work/production/Broome_2008/build_boundary_27255.py
===================================================================
--- /anuga_work/production/Broome_2008/build_boundary_27255.py	(revision 7159)
+++ 	(revision )
@@ -1,146 +1,0 @@
-"""
-Script for building boundary to run tsunami inundation scenario for Broome,
-WA, Australia.  The boundary is based on the National Hazard Map.
-
-Input: order_filename from project.py
-Output: creates a sts file and csv files stored in project.boundaries_dir_event.
-The run_broome.py is reliant on the output of this script.
-"""
-import os
-from os import sep
-import project
-
-from anuga.utilities.numerical_tools import ensure_numeric
-from Scientific.IO.NetCDF import NetCDFFile
-from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
-    compress,zeros,fabs,allclose,ones
-from anuga.utilities.polygon import inside_polygon,read_polygon
-from time import localtime, strftime, gmtime
-from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
-
-try:
-    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
-except:
-    print 'Cannot import pylab plotting will not work. Csv files are still created'
-
-#--------------------------------------------------------------------------
-# Create sts boundary from mux2files 
-#--------------------------------------------------------------------------
-
-dir=os.path.join(project.muxhome,'mux')
-
-# Refer to event_027255.list in home+state+sep+event08 for event details
-# taken from David's event list for 27255
-urs_filenames = [
-    os.path.join(dir,'Java-0017-z.grd'),  
-    os.path.join(dir,'Java-0018-z.grd'),
-    os.path.join(dir,'Java-0019-z.grd'),
-    os.path.join(dir,'Java-0022-z.grd'),  
-    os.path.join(dir,'Java-0023-z.grd'),  
-    os.path.join(dir,'Java-0024-z.grd'),  
-    os.path.join(dir,'Java-0027-z.grd'),  
-    os.path.join(dir,'Java-0028-z.grd'),  
-    os.path.join(dir,'Java-0031-z.grd'),  
-    os.path.join(dir,'Java-0032-z.grd')]  
- 
-print 'number of sources', len(urs_filenames)
-
-if len(urs_filenames) == 10: #change per event
-
-    weight_factor = 30.84050 #change per event
-    weights=weight_factor*ones(len(urs_filenames),Float)
-
-    scenario_name=project.scenario_name
-    order_filename=os.path.join(project.order_filename_dir)
-
-    print 'reading', order_filename
-    # Create ordered sts file
-    print 'creating sts file'
-
-    urs2sts(urs_filenames,
-            basename_out=os.path.join(project.boundaries_dir_event,scenario_name),
-            ordering_filename=order_filename,
-            weights=weights,
-            mean_stage=project.tide,
-            verbose=True)
-else:
-    print 'number of sources do not match event.list'
-
-#------------------------------------------------------------------------------
-# Get gauges (timeseries of index points)
-#------------------------------------------------------------------------------
-def get_sts_gauge_data(filename,verbose=False):
-    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
-        compress,zeros,fabs,take
-    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=transpose(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 Maxium wave height throughout timeseries at each index point
-#------------------------------------------------------------------------------
-
-    maxname = 'max_sts_stage.csv'
-    fid_max = open(project.boundaries_dir_event+sep+maxname,'w')
-    s = 'index, x, y, max_stage \n'
-    fid_max.write(s)    
-    for j in range(len(x)):
-        index= permutation[j]
-        stage = quantities['stage'][:,j]
-        xmomentum = quantities['xmomentum'][:,j]
-        ymomentum = quantities['ymomentum'][:,j]
-
-        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
-        fid_max.write(s)
-      
-#------------------------------------------------------------------------------
-# Get Minium wave height throughout timeseries at each index point
-#------------------------------------------------------------------------------
-
-    minname = 'min_sts_stage.csv'
-    fid_min = open(project.boundaries_dir_event+sep+minname,'w')
-    s = 'index, x, y, max_stage \n'
-    fid_min.write(s)    
-    for j in range(len(x)):
-        index= permutation[j]
-        stage = quantities['stage'][:,j]
-        xmomentum = quantities['xmomentum'][:,j]
-        ymomentum = quantities['ymomentum'][:,j]
-
-        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))
-        fid_min.write(s)
-        
-
-
-        fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w')
-        s = 'time, stage, xmomentum, ymomentum \n'
-        fid_sts.write(s)
-
-#------------------------------------------------------------------------------
-# End of the get gauges
-#------------------------------------------------------------------------------
-        for k in range(len(time)-1):
-            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
-            fid_sts.write(s)
-
-        fid_sts.close()      
-
-    fid.close()
-
-
-    return quantities,elevation,time
-
-quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False)
-
-print len(elevation), len(quantities['stage'][0,:])
-
Index: uga_work/production/Broome_2008/build_boundary_27283.py
===================================================================
--- /anuga_work/production/Broome_2008/build_boundary_27283.py	(revision 7159)
+++ 	(revision )
@@ -1,185 +1,0 @@
-"""
-Script for building boundary to run tsunami inundation scenario for Broome,
-WA, Australia.  The boundary is based on the National Hazard Map.
-
-Input: order_filename from project.py
-Output: creates a sts file and csv files stored in project.boundaries_dir_event.
-The run_broome.py is reliant on the output of this script.
-"""
-import os
-from os import sep
-import project
-
-from anuga.utilities.numerical_tools import ensure_numeric
-from Scientific.IO.NetCDF import NetCDFFile
-from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
-    compress,zeros,fabs,allclose,ones
-from anuga.utilities.polygon import inside_polygon,read_polygon
-from time import localtime, strftime, gmtime
-from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
-
-try:
-    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
-except:
-    print 'Cannot import pylab plotting will not work. Csv files are still created'
-
-#--------------------------------------------------------------------------
-# Create sts boundary from mux2files 
-#--------------------------------------------------------------------------
-
-dir=os.path.join(project.muxhome,'mux')
-
-# Refer to event_027255.list in home+state+sep+event08 for event details
-# taken from David's event list for 27255
-urs_filenames = [
-    os.path.join(dir,'Java-0000-z.grd'),
-    os.path.join(dir,'Java-0001-z.grd'),
-    os.path.join(dir,'Java-0002-z.grd'),
-    os.path.join(dir,'Java-0003-z.grd'),
-    os.path.join(dir,'Java-0005-z.grd'), 
-    os.path.join(dir,'Java-0006-z.grd'), 
-    os.path.join(dir,'Java-0007-z.grd'), 
-    os.path.join(dir,'Java-0008-z.grd'), 
-    os.path.join(dir,'Java-0010-z.grd'), 
-    os.path.join(dir,'Java-0011-z.grd'), 
-    os.path.join(dir,'Java-0012-z.grd'), 
-    os.path.join(dir,'Java-0013-z.grd'),  
-    os.path.join(dir,'Java-0015-z.grd'),  
-    os.path.join(dir,'Java-0016-z.grd'),  
-    os.path.join(dir,'Java-0017-z.grd'),  
-    os.path.join(dir,'Java-0018-z.grd'),  
-    os.path.join(dir,'Java-0020-z.grd'),  
-    os.path.join(dir,'Java-0021-z.grd'),  
-    os.path.join(dir,'Java-0022-z.grd'),  
-    os.path.join(dir,'Java-0023-z.grd'),  
-    os.path.join(dir,'Java-0025-z.grd'),  
-    os.path.join(dir,'Java-0026-z.grd'),  
-    os.path.join(dir,'Java-0027-z.grd'),  
-    os.path.join(dir,'Java-0028-z.grd'),  
-    os.path.join(dir,'Java-0029-z.grd'),  
-    os.path.join(dir,'Java-0030-z.grd'),  
-    os.path.join(dir,'Java-0031-z.grd'),  
-    os.path.join(dir,'Java-0032-z.grd'),  
-    os.path.join(dir,'Java-0033-z.grd'),  
-    os.path.join(dir,'Java-0034-z.grd'),  
-    os.path.join(dir,'Java-0035-z.grd'),  
-    os.path.join(dir,'Java-0036-z.grd'),  
-    os.path.join(dir,'Java-0037-z.grd'),  
-    os.path.join(dir,'Java-0038-z.grd'),  
-    os.path.join(dir,'Java-0039-z.grd'),  
-    os.path.join(dir,'Java-0040-z.grd'),  
-    os.path.join(dir,'Java-0041-z.grd'),  
-    os.path.join(dir,'Java-0042-z.grd'),  
-    os.path.join(dir,'Java-0043-z.grd'),  
-    os.path.join(dir,'Java-0044-z.grd'),  
-    os.path.join(dir,'Java-0045-z.grd'),  
-    os.path.join(dir,'Java-0046-z.grd'),  
-    os.path.join(dir,'Java-0047-z.grd'),
-    os.path.join(dir,'Java-0048-z.grd')]    
- 
-print 'number of sources', len(urs_filenames)
-
-if len(urs_filenames) == 44: #change per event
-
-# AS per David Burbidge email on Friday 4th July the mag 9.3 event
-# has 1m worth of slip on each sub fault therefore mutliple each unit
-# source by the slip (10.4544) and sum the 44 time series together to
-# get the time series for this event at the points on your boundary. 
-
-    weight_factor = 10.4544 #change per event
-    weights=weight_factor*ones(len(urs_filenames),Float)
-
-    scenario_name=project.scenario_name
-    order_filename=os.path.join(project.order_filename_dir)
-
-    print 'reading', order_filename
-    # Create ordered sts file
-    print 'creating sts file'
-
-    urs2sts(urs_filenames,
-            basename_out=os.path.join(project.boundaries_dir_event,scenario_name),
-            ordering_filename=order_filename,
-            weights=weights,
-            mean_stage=project.tide,
-            verbose=True)
-else:
-    print 'number of sources do not match event.list'
-
-#------------------------------------------------------------------------------
-# Get gauges (timeseries of index points)
-#------------------------------------------------------------------------------
-def get_sts_gauge_data(filename,verbose=False):
-    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
-        compress,zeros,fabs,take
-    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=transpose(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 Maxium wave height throughout timeseries at each index point
-#------------------------------------------------------------------------------
-
-    maxname = 'max_sts_stage.csv'
-    fid_max = open(project.boundaries_dir_event+sep+maxname,'w')
-    s = 'index, x, y, max_stage \n'
-    fid_max.write(s)    
-    for j in range(len(x)):
-        index= permutation[j]
-        stage = quantities['stage'][:,j]
-        xmomentum = quantities['xmomentum'][:,j]
-        ymomentum = quantities['ymomentum'][:,j]
-
-        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
-        fid_max.write(s)
-      
-#------------------------------------------------------------------------------
-# Get Minium wave height throughout timeseries at each index point
-#------------------------------------------------------------------------------
-
-    minname = 'min_sts_stage.csv'
-    fid_min = open(project.boundaries_dir_event+sep+minname,'w')
-    s = 'index, x, y, max_stage \n'
-    fid_min.write(s)    
-    for j in range(len(x)):
-        index= permutation[j]
-        stage = quantities['stage'][:,j]
-        xmomentum = quantities['xmomentum'][:,j]
-        ymomentum = quantities['ymomentum'][:,j]
-
-        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))
-        fid_min.write(s)
-        
-
-
-        fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w')
-        s = 'time, stage, xmomentum, ymomentum \n'
-        fid_sts.write(s)
-
-#------------------------------------------------------------------------------
-# End of the get gauges
-#------------------------------------------------------------------------------
-        for k in range(len(time)-1):
-            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
-            fid_sts.write(s)
-
-        fid_sts.close()      
-
-    fid.close()
-
-
-    return quantities,elevation,time
-
-quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False)
-
-print len(elevation), len(quantities['stage'][0,:])
-
Index: uga_work/production/Broome_2008/build_broome.py
===================================================================
--- /anuga_work/production/Broome_2008/build_broome.py	(revision 7159)
+++ 	(revision )
@@ -1,145 +1,0 @@
-"""
-Script for building the elevation data to run tsunami inundation scenario
-for Broome, WA, Australia.
-
-Input: elevation data from project.py
-Output: pts file stored in project.topographies_dir 
-The run_broome.py is reliant on the output of this script.
-
-Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
-"""
-
-#------------------------------------------------------------------------------
-# Import necessary modules
-#------------------------------------------------------------------------------
-
-# Standard modules
-from os import sep
-from os.path import dirname, basename
-from os import mkdir, access, F_OK
-from shutil import copy
-import time
-import sys
-
-# Related major packages
-from anuga.shallow_water import Domain
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.geospatial_data.geospatial_data import *
-from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
-
-# Application specific imports
-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_time_dir,__file__, 
-               dirname(project.__file__)+sep+ project.__name__+'.py' )
-
-start_screen_catcher(project.output_build_time_dir)
-
-print 'USER:    ', project.user
-
-#-------------------------------------------------------------------------------
-# Preparation of topographic data
-#
-# Convert ASC 2 DEM 2 PTS using source data and store result in source data
-# Do for coarse and fine data
-# Fine pts file to be clipped to area of interest
-#-------------------------------------------------------------------------------
-print"project.poly_all",project.poly_all
-print"project.combined_dir_name",project.combined_dir_name
-
-# input elevation directory filenames
-onshore_in_dir_name = project.onshore_in_dir_name
-
-coast_in_dir_name = project.coast_in_dir_name
-
-offshore_in_dir_name = project.offshore_in_dir_name
-offshore_in_dir_name1 = project.offshore_in_dir_name1
-offshore_in_dir_name2 = project.offshore_in_dir_name2
-offshore_in_dir_name3 = project.offshore_in_dir_name3
-offshore_in_dir_name4 = project.offshore_in_dir_name4
-offshore_in_dir_name5 = project.offshore_in_dir_name5
-offshore_in_dir_name6 = project.offshore_in_dir_name6
-
-# output elevation directory filenames
-onshore_dir_name = project.onshore_dir_name
-
-coast_dir_name = project.coast_dir_name
-
-offshore_dir_name = project.offshore_dir_name
-offshore_dir_name1 = project.offshore_dir_name1
-offshore_dir_name2 = project.offshore_dir_name2
-offshore_dir_name3 = project.offshore_dir_name3
-offshore_dir_name4 = project.offshore_dir_name4
-offshore_dir_name5 = project.offshore_dir_name5
-offshore_dir_name6 = project.offshore_dir_name6
-
-# creates DEM from asc data
-print "creates DEMs from asc data"
-convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
-
-convert_dem_from_ascii2netcdf(offshore_in_dir_name1, basename_out=offshore_dir_name1, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(offshore_in_dir_name2, basename_out=offshore_dir_name2, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(offshore_in_dir_name3, basename_out=offshore_dir_name3, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(offshore_in_dir_name4, basename_out=offshore_dir_name4, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)
-convert_dem_from_ascii2netcdf(offshore_in_dir_name6, basename_out=offshore_dir_name6, use_cache=True, verbose=True)
-
-#creates pts file for onshore DEM
-print "creates pts file for onshore DEM"
-dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name1 ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name2 ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name3 ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name4 ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
-dem2pts(offshore_dir_name6 ,use_cache=True,verbose=True)
-
-###creates pts file for island DEM
-##dem2pts(island_dir_name, use_cache=True, verbose=True)
-##dem2pts(island_dir_name1, use_cache=True, verbose=True)
-##dem2pts(island_dir_name2, use_cache=True, verbose=True)
-##dem2pts(island_dir_name3, use_cache=True, verbose=True)
-##dem2pts(island_dir_name4, use_cache=True, verbose=True)
-##dem2pts(island_dir_name5, use_cache=True, verbose=True)
-##dem2pts(island_dir_name6, use_cache=True, verbose=True)
-
-#-------------------------------------------------------------------------------
-# Combine datasets into project.combined_dir_name
-#-------------------------------------------------------------------------------
-
-print'create Geospatial data1 objects from topographies'
-G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
-print'create Geospatial data6 objects from topographies'
-G_coast = Geospatial_data(file_name = coast_in_dir_name)
-print'create Geospatial data7 objects from topographies'
-G_off = Geospatial_data(file_name = offshore_in_dir_name)
-print'create Geospatial data8 objects from topographies'
-G_off1 = Geospatial_data(file_name = offshore_dir_name1 + '.pts')
-print'create Geospatial data9 objects from topographies'
-G_off2 = Geospatial_data(file_name = offshore_dir_name2 + '.pts')
-print'create Geospatial data10 objects from topographies'
-G_off3 = Geospatial_data(file_name = offshore_dir_name3 + '.pts')
-print'create Geospatial data10 objects from topographies'
-G_off4 = Geospatial_data(file_name = offshore_dir_name4 + '.pts')
-print'create Geospatial data10 objects from topographies'
-G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
-print'create Geospatial data10 objects from topographies'
-G_off6 = Geospatial_data(file_name = offshore_dir_name6 + '.pts')
-print'add all geospatial objects'
-G = G1 + G_coast + G_off + G_off1 + G_off2 + G_off3  + G_off4 + G_off5 + G_off6
-print'clip combined geospatial object by bounding polygon'
-G_clipped = G.clip(project.poly_all)
-
-print'export combined DEM file'
-if access(project.topographies_dir,F_OK) == 0:
-    mkdir (project.topographies_dir)
-G_clipped.export_points_file(project.combined_dir_name + '.pts')
-#G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
-
-import sys
-sys.exit()
Index: uga_work/production/Broome_2008/export_results_all.py
===================================================================
--- /anuga_work/production/Broome_2008/export_results_all.py	(revision 7159)
+++ 	(revision )
@@ -1,121 +1,0 @@
-"""
-Generates ascii grids of nominated areas -
-Input: sww file from run_perth.py
-       boundaries for grids from project.py
-Outputs: ascii grids of specified variables
-Stored in the 'outputs_dir' folder for respective .sww file
-
-Note:
-If producing a grid for the enitre extent cellsize should be greater than 30m
-If producing grids for inundation area resolution should be greater than mesh (ie ~22m)
-"""
-
-import project, os
-import sys
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.shallow_water.data_manager import sww2dem
-from os import sep
-
-
-directory = project.output_dir
-
-#time_dir = '20080526_104946_run_final_0.6_test_kvanputt'
-#time_dir = '20080530_170833_run_final_0.6_exmouth_kvanputt'
-time_dir1 = '20081002_171412_run_final_0.6_27255_alpha0.1_rwilson'
-time_dirs = [time_dir1]
-
-#cellsize = 20
-cellsize = 30
-
-#timestep = None    # None means no timestep!
-timestep = 0
-
-######
-# Set the special areas of interest.  If none, do: area='All'
-######
-#area = ['Geordie', 'Sorrento', 'Fremantle', 'Rockingham']  # strings must match keys in var_equations below
-area = ['All']      # 'All' means no special areas - the whole thing
-
-######
-# Define allowed variable names and associated equations to generate values.
-# This would not normally change.
-######
-var_equations = {'stage':     'stage',
-                 'momentum':  '(xmomentum**2 + ymomentum**2)**0.5',
-                 'depth':     'stage-elevation',
-                 'speed':     '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)',
-                 'elevation': 'elevation' }
-
-# one or more key strings from var_equations above
-var = ['stage', 'elevation']
-
-######
-# Start running the various conversions we require.
-######
-for time_dir in time_dirs:
-    name1 = directory+time_dir+sep+project.scenario_name
-#    name2 = directory+time_dir+sep+'sww2'+sep+project.scenario_name+'_time_39900_0' #need to get assistance on how to make this into anything
-    names = [name1]
-    
-    for name in names:
-        for which_area in area:
-            if which_area == 'All':
-                easting_min = None
-                easting_max = None
-                northing_min = None
-                northing_max = None
-            else:
-                try:
-                    easting_min = eval('project.xmin%s' % which_area)
-                    easting_max = eval('project.xmax%s' % which_area)
-                    northing_min = eval('project.ymin%s' % which_area)
-                    northing_max = eval('project.ymax%s' % which_area)
-                except AttributeError:
-                    print 'Unrecognized area name: %s' % which_area
-                    break
-### This is replaced by the single 'else' case above
-##            elif which_area == 'Geordie':
-##                easting_min = project.xminGeordie
-##                easting_max = project.xmaxGeordie
-##                northing_min = project.yminGeordie
-##                northing_max = project.ymaxGeordie
-##            elif which_area == 'Sorrento':
-##                easting_min = project.xminSorrento
-##                easting_max = project.xmaxSorrento
-##                northing_min = project.yminSorrento
-##                northing_max = project.ymaxSorrento
-##            elif which_area == 'Fremantle':
-##                easting_min = project.xminFremantle
-##                easting_max = project.xmaxFremantle
-##                northing_min = project.yminFremantle
-##                northing_max = project.ymaxFremantle
-##            elif which_area == 'Rockingham':
-##                easting_min = project.xminRockingham
-##                easting_max = project.xmaxRockingham
-##                northing_min = project.yminRockingham
-##                northing_max = project.ymaxRockingham
-##            else:
-##                print 'Unrecognized area name: %s' % which_area
-##                break
-              
-            for which_var in var:
-                if which_var not in var_equations:
-                    print 'Unrecognized variable name: %s' % which_var
-                    break
-
-                outname = name + which_area + '_' + which_var
-                quantityname = var_equations[which_var]
-
-                print 'start sww2dem: time_dir=%s' % time_dir
-                
-                sww2dem(name, basename_out = outname,
-                            quantity = quantityname,
-                            timestep = timestep,
-                            cellsize = cellsize,      
-                            easting_min = easting_min,
-                            easting_max = easting_max,
-                            northing_min = northing_min,
-                            northing_max = northing_max,        
-                            reduction = max, 
-                            verbose = True,
-                            format = 'asc')
Index: uga_work/production/Broome_2008/export_results_max.py
===================================================================
--- /anuga_work/production/Broome_2008/export_results_max.py	(revision 7159)
+++ 	(revision )
@@ -1,115 +1,0 @@
-"""
-Generates ascii grids of nominated areas -
-Input: sww file from run_broome.py
-       boundaries for grids from project.py
-Outputs: ascii grids of specified variables
-Stored in the 'outputs_dir' folder for respective .sww file
-
-Note:
-If producing a grid for the enitre extent cellsize should be greater than 30m
-If producing grids for inundation area resolution should be greater than mesh (ie ~22m)
-"""
-
-import project, os
-import sys
-from anuga.lib.maxasc.maxasc import MaxAsc
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.shallow_water.data_manager import sww2dem
-from os import sep
-
-
-directory = project.output_dir
-
-
-time_dir = '20090106_151505_run_final_0.6_27283_alpha0.1_jgriffin'
-
-time_dirs = [time_dir]# , time_dir3, time_dir4, time_dir5, time_dir6]
-  
-
-#cellsize = 20
-cellsize = 250
-
-timestep = None    # None means no timestep!
-#timestep = 0
-
-######
-# Set the special areas of interest.  If none, do: area='All'
-######
-
-#area = ['Bunbury', 'Busselton']  # strings must match keys in var_equations below
-area = ['All']      # 'All' means no special areas - the whole thing
-
-######
-# Define allowed variable names and associated equations to generate values.
-# This would not normally change.
-######
-var_equations = {'stage':     'stage',
-                 'momentum':  '(xmomentum**2 + ymomentum**2)**0.5',
-                 'depth':     'stage-elevation',
-                 'speed':     '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)',
-                 'elevation': 'elevation' }
-
-# one or more key strings from var_equations above
-var = ['stage','depth','speed']
-
-######
-# Start running the various conversions we require.
-######
-
-for which_var in var:
-    if which_var not in var_equations:
-        print 'Unrecognized variable name: %s' % which_var
-        break
-
-    for which_area in area:
-        if which_area == 'All':
-            easting_min = None
-            easting_max = None
-            northing_min = None
-            northing_max = None
-        else:
-            try:
-                easting_min = eval('project.xmin%s' % which_area)
-                easting_max = eval('project.xmax%s' % which_area)
-                northing_min = eval('project.ymin%s' % which_area)
-                northing_max = eval('project.ymax%s' % which_area)
-            except AttributeError:
-                print 'Unrecognized area name: %s' % which_area
-                break
-
-        for time_dir in time_dirs:
-
-            name1 = directory+time_dir+sep+project.scenario_name
-            name2 = directory+time_dir+sep+project.scenario_name+'_time_63900_0'
-          
-            names = [name1, name2]
-      
-            asc_name = []   
-               
-            for name in names:
-                         
-                outname = name + '_' + which_area + '_' + which_var
-                quantityname = var_equations[which_var]
-
-                print 'start sww2dem: time_dir=%s' % time_dir
-                
-                sww2dem(name, basename_out = outname,
-                            quantity = quantityname,
-                            timestep = timestep,
-                            cellsize = cellsize,      
-                            easting_min = easting_min,
-                            easting_max = easting_max,
-                            northing_min = northing_min,
-                            northing_max = northing_max,        
-                            reduction = max, 
-                            verbose = True,
-                            format = 'asc')
-
-                asc_name.append(outname + '.asc')
-    
-            maxasc_outname = directory+time_dir+sep+project.scenario_name+'_'+which_area+'_'+which_var+'_max.asc'
-
-            print 'max asc outname ', maxasc_outname
-            print 'asc_name ', str(asc_name)
-            
-            MaxAsc(maxasc_outname, asc_name)
Index: uga_work/production/Broome_2008/get_timeseries.py
===================================================================
--- /anuga_work/production/Broome_2008/get_timeseries.py	(revision 7159)
+++ 	(revision )
@@ -1,40 +1,0 @@
-"""
-Generate time series of nominated "gauges".
-Input: sww file from run_perth.py
-       gauges project.gauge_filename
-Output: csv files stage,speed,depth,elevation over time
-Stored in the 'outputs_dir' folder for respective .sww file
-
-Note:
-Can run different sww files at the same time
-make sure if there is a second sww file that it is placed into a folder called sww2
-Can run different gauges at the same time - ie testing boundary index point
-"""
-
-from os import sep
-from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
-import project
-
-timestamp1='20080815_103442_run_final_0.0_polyline_alpha0.1_kvanputt'
-timestamp2='20080815_103336_run_final_0.6_polyline_alpha0.1_kvanputt'
-
-timestamps = [timestamp1, timestamp2]
-for timestamp in timestamps:
-
-    filename1=project.output_dir+timestamp+sep+project.scenario_name+'.sww'
-    filename2=project.output_dir+timestamp+sep+'sww2'+sep+project.scenario_name+'_time_39900_0.sww'
-
-    filenames = [filename1, filename2]
-    for filename in filenames:
-
-        gauges = [project.gauges_dir_name, project.gauges_dir_name2]
-        for gauge in gauges:
-            
-            sww2csv_gauges(filename,gauge,
-                            quantities = ['stage','speed','depth','elevation'],
-                            verbose=True)
-                        
-
-    
-
-
Index: uga_work/production/Broome_2008/locate_indigeneous.py
===================================================================
--- /anuga_work/production/Broome_2008/locate_indigeneous.py	(revision 7159)
+++ 	(revision )
@@ -1,41 +1,0 @@
-#
-# script to determine which communities are located in study region
-# 31 October 2006
-#
-
-def community_in_extent(filename,extent,zone,fid_out):
-    fid = open(filename)
-    lines = fid.readlines()
-    fid.close()
-
-    line1 = lines[0]
-    line11 = line1.split(',')
-    for i in range(len(line11)):
-        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
-        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
-        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
-
-    from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
-    from anuga.utilities.polygon import plot_polygons, is_inside_polygon
-    
-    for line in lines[1:]:
-        fields = line.split(',')
-        point, zone1 = convert_from_latlon_to_utm(latitudes=[fields[lat_index]],\
-                                                 longitudes=[fields[lon_index]])
-        if zone1 == zone:
-            v = is_inside_polygon(point,extent,verbose=False)
-            if v == True:
-                p = point[0]
-                s = '%s, %.2f, %.2f\n' %(fields[name_index], p[0], p[1])
-                fid_out.write(s)
-    return
-
-import project
-
-filename = project.community_filename
-filename_out = project.community_perth
-fid_out = open(filename_out, 'w')
-extent = project.polyAll
-zone = 50
-
-community_in_extent(filename,extent,zone,fid_out)
Index: uga_work/production/Broome_2008/plot_timeseries.py
===================================================================
--- /anuga_work/production/Broome_2008/plot_timeseries.py	(revision 7159)
+++ 	(revision )
@@ -1,32 +1,0 @@
-
-
-from os import sep
-import project
-from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
-
-import pylab
-
-timestamp1='20080815_103336_run_final_0.6_polyline_alpha0.1_kvanputt'
-timestamp2='20080815_103442_run_final_0.0_polyline_alpha0.1_kvanputt'
-
-timestamps = [timestamp1, timestamp2]
-for timestamp in timestamps:
-
-    input_dir_name1=project.output_dir+timestamp+sep
-    input_dir_name2=project.output_dir+timestamp+sep+'perth_time_39900'+sep
-
-    input_dir_names = [input_dir_name1, input_dir_name2]
-    for input_dir_name in input_dir_names:
-
-        print input_dir_name
-
-        csv2timeseries_graphs(directories_dic={input_dir_name:['test',0, 0]},
-                                output_dir=input_dir_name,
-                                base_name='gauge_',
-                                plot_numbers='',
-                                quantities=['stage','speed','depth'],
-                                extra_plot_name='',
-                                assess_all_csv_files=True,                            
-                                create_latex=False,
-                                verbose=True)
-
Index: uga_work/production/Broome_2008/project.py
===================================================================
--- /anuga_work/production/Broome_2008/project.py	(revision 7159)
+++ 	(revision )
@@ -1,248 +1,0 @@
-"""Common filenames and locations for elevation, meshes and outputs.
-This script is the heart of all scripts in folder
-"""
-#------------------------------------------------------------------------------
-# Import necessary modules
-#------------------------------------------------------------------------------
-
-from os import sep, environ, getenv, getcwd
-from os.path import expanduser
-import sys
-from time import localtime, strftime, gmtime
-from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
-from anuga.utilities.system_tools import get_user_name, get_host_name
-from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
-
-#------------------------------------------------------------------------------
-# Directory setup
-#------------------------------------------------------------------------------
-##Note: INUNDATIONHOME is the inundation directory, not the data directory.
-
-home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
-muxhome = getenv('MUXHOME')
-user = get_user_name()
-host = get_host_name()
-
-##time stuff
-time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
-gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
-build_time = time+'_build'
-run_time = time+'_run'
-print 'gtime: ', gtime
-
-#------------------------------------------------------------------------------
-# Initial Conditions
-#------------------------------------------------------------------------------
-
-##Changed to reflect the modeled community (also sets up the directory system)
-state = 'western_australia'
-scenario_name = 'broome'
-scenario = 'broome_tsunami_scenario_2008'
-##One or all can be changed each time the run_scenario script is excuted
-tide = 0.6
-event_number = 27255
-#event_number = 27283
-alpha = 0.1
-friction=0.01
-starttime=0
-finaltime=80000 
-
-setup='final'  ## Can replace with trial or basic, this action will thin the mesh
-               ## so that the run script will take less time (hence less acurate) 
-
-if setup =='trial':
-    print'trial'
-    res_factor=10
-    time_thinning=48
-    yieldstep=240
-if setup =='basic': 
-    print'basic'
-    res_factor=4
-    time_thinning=12
-    yieldstep=120
-if setup =='final': 
-    print'final'
-    res_factor=1
-    time_thinning=4
-    yieldstep=60
-
-#------------------------------------------------------------------------------
-# Output Filename
-#------------------------------------------------------------------------------
-##Important to distiquish each run
-dir_comment='_'+setup+'_'+str(tide)+'_'+str(event_number)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
-
-#------------------------------------------------------------------------------
-# INPUT DATA
-#------------------------------------------------------------------------------
-
-##ELEVATION DATA## - used in build_broome.py
-## onshore data: format ascii grid with accompaning projection file
-onshore_name = 'town_topo_10m' 
-## island: format ascii grid with accompaning projection file
-#island_name = 'rott_dli_ext' 
-#island_name1 = 'gard_dli_ext'
-#island_name2 = 'carnac_island_dli_ext'
-#island_name3 = 'penguin_dli_ext'
-## coastline: format x,y,elevation (with title)
-coast_name = 'Broome_coastline.txt'
-## bathymetry: format x,y,elevation (with title)
-offshore_name = 'Broome_Bathymetry.txt'
-offshore_name1 = 'cable_250m'
-offshore_name2 = 'inferred_north'
-offshore_name3 = 'inferred_south'
-offshore_name4 = 'north_250m'
-offshore_name5 = 'other_topo_250m'
-offshore_name6 = 'south_250m'
-
-##GAUGES## - used in get_timeseries.py
-#gauge_name = 'Broome.csv'
-#gauge_name2 = 'thinned_MGA50.csv'
-
-##BOUNDING POLYGON## -used in build_boundary.py and run_broome.py respectively
-#NOTE: when files are put together must be in sequence - for ease go clockwise!
-#Check the run_broome.py for boundary_tags
-##thinned ordering file from Hazard Map: format index,latitude,longitude (with title)
-order_filename = 'thinned_boundary_ordering.csv'
-##landward bounding points
-landward = 'landward_bounding_polygon.csv'
-
-#------------------------------------------------------------------------------
-# OUTPUT ELEVATION DATA
-#------------------------------------------------------------------------------
-##Output filename for elevation
-## its a combination of all the data put together (utilisied in build_boundary)
-combined_name ='broome_combined_elevation'
-combined_smaller_name = 'broome_combined_elevation_smaller'
-
-#------------------------------------------------------------------------------
-# Directory Structure
-#------------------------------------------------------------------------------
-anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
-topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
-topographies_dir = anuga_dir+'topographies'+sep
-polygons_dir = anuga_dir+'polygons'+sep
-tide_dir = anuga_dir+'tide_data'+sep
-boundaries_dir = anuga_dir+'boundaries'+ sep 
-output_dir = anuga_dir+'outputs'+sep
-gauges_dir = anuga_dir+'gauges'+sep
-meshes_dir = anuga_dir+'meshes'+sep
-
-#------------------------------------------------------------------------------
-# Location of input and output Data
-#------------------------------------------------------------------------------
-##where the input data sits
-onshore_in_dir_name = topographies_in_dir + onshore_name
-coast_in_dir_name = topographies_in_dir + coast_name
-offshore_in_dir_name = topographies_in_dir + offshore_name
-offshore_in_dir_name1 = topographies_in_dir + offshore_name1
-offshore_in_dir_name2 = topographies_in_dir + offshore_name2
-offshore_in_dir_name3 = topographies_in_dir + offshore_name3
-offshore_in_dir_name4 = topographies_in_dir + offshore_name4
-offshore_in_dir_name5 = topographies_in_dir + offshore_name5
-offshore_in_dir_name6 = topographies_in_dir + offshore_name6
-
-##where the output data sits
-onshore_dir_name = topographies_dir + onshore_name
-coast_dir_name = topographies_dir + coast_name
-offshore_dir_name = topographies_dir + offshore_name
-offshore_dir_name1 = topographies_dir + offshore_name1
-offshore_dir_name2 = topographies_dir + offshore_name2
-offshore_dir_name3 = topographies_dir + offshore_name3
-offshore_dir_name4 = topographies_dir + offshore_name4
-offshore_dir_name5 = topographies_dir + offshore_name5
-offshore_dir_name6 = topographies_dir + offshore_name6
-
-##where the combined file sits
-combined_dir_name = topographies_dir + combined_name
-combined_smaller_name_dir = topographies_dir + combined_smaller_name
-
-##where the mesh sits (this is created during the run_broome.py)
-meshes_dir_name = meshes_dir + scenario_name+'.msh'
-
-#where the boundary order files sit (this is used within build_boundary.py)
-order_filename_dir = boundaries_dir + order_filename
-
-#where the landward points of boundary extent sit (this is used within run_broome.py)
-landward_dir = boundaries_dir + landward
-
-#where the event sts files sits (this is created during the build_boundary.py)
-boundaries_dir_event = boundaries_dir + str(event_number) + sep
-boundaries_dir_mux = muxhome
-
-#where the directory of the output filename sits
-output_build_time_dir = output_dir+build_time+dir_comment+sep #used for build_broome.py
-output_run_time_dir = output_dir+run_time+dir_comment+sep #used for run_broome.py
-output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
-
-#where the directory of the gauges sit
-#gauges_dir_name = gauges_dir + gauge_name #used for get_timeseries.py
-#gauges_dir_name2 = gauges_dir + gauge_name2 #used for get_timeseries.py
-
-#------------------------------------------------------------------------------
-# Interior region definitions
-#------------------------------------------------------------------------------
-
-#Land, to set the stage/water to be offcoast only
-poly_mainland = read_polygon(polygons_dir+'initial_conditions_mainland.csv')
-
-#Initial bounding polygon for data clipping 
-poly_all = read_polygon(polygons_dir+'poly_all.csv')
-res_poly_all = 100000*res_factor
-
-#Area of Interest 1 (Broome)
-poly_aoi1 = read_polygon(polygons_dir+'Area_Of_Interest.csv')
-res_aoi1 = 500*res_factor
-
-#Area of Significance 1 (Broome)
-poly_aos1 = read_polygon(polygons_dir+'Area_Of_Significance.csv')
-res_aos1 = 1000*res_factor
-
-#Shallow water 1
-poly_sw1 = read_polygon(polygons_dir+'Shallow_Water.csv')
-res_sw1 = 25000*res_factor
-
-##Deep water (land of Rottnest, ANUGA does not do donuts!)
-#poly_dw1 = read_polygon(polygons_dir+'rottnest_internal.csv')
-#res_dw1 = 100000*res_factor
-
-# Combined all regions, must check that all are included!
-interior_regions = [[poly_aoi1,res_aoi1],
-                    [poly_aos1,res_aos1],
-                    [poly_sw1,res_sw1]]
-
-    
-trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
-print 'min estimated number of triangles', trigs_min
-   
-
-###------------------------------------------------------------------------------
-### Clipping regions for export to asc and regions for clipping data
-###------------------------------------------------------------------------------
-##
-###Geordie Bay extract ascii grid
-##xminGeordie = 358000
-##xmaxGeordie = 362000
-##yminGeordie = 6458500
-##ymaxGeordie = 6461000
-##
-###Sorrento extract ascii grid
-##xminSorrento = 379000
-##xmaxSorrento = 382500
-##yminSorrento = 6477000
-##ymaxSorrento = 6480000
-##
-###Fremantle extract ascii grid
-##xminFremantle = 376000
-##xmaxFremantle = 388000
-##yminFremantle = 6449000
-##ymaxFremantle = 6461000
-##
-###Rockingham extract ascii grid
-##xminRockingham = 373500
-##xmaxRockingham = 385500
-##yminRockingham = 6424000
-##ymaxRockingham = 6433000
-##
-##
Index: uga_work/production/Broome_2008/run_broome.py
===================================================================
--- /anuga_work/production/Broome_2008/run_broome.py	(revision 7159)
+++ 	(revision )
@@ -1,250 +1,0 @@
-"""Script for running tsunami inundation scenario for Broome, WA, Australia.
-
-The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data is combiled into a pts file through build_broome.py
-and a simulated tsunami is generated through an sts file from build_boundary.py
-
-Input: sts file (build_boundary.py for resepective event)
-       pts file (build_broome.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
-from os import sep
-import os
-from os.path import dirname, basename
-from os import mkdir, access, F_OK
-from shutil import copy
-import time
-import sys
-
-# Related major packages
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Dirichlet_boundary
-from anuga.shallow_water import File_boundary
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Field_boundary
-from Numeric import allclose
-from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
-from anuga.pmesh.mesh_interface import create_mesh_from_regions
-from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
-from anuga_parallel.parallel_abstraction import get_processor_name
-from anuga.caching import myhash
-from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
-from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
-from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
-from polygon import Polygon_function
-    
-# Application specific imports
-import project  # Definition of file names and polygons
-numprocs = 1
-myid = 0
-
-def run_model(**kwargs):
-    
-    #------------------------------------------------------------------------------
-    # Copy scripts to time stamped output directory and capture screen
-    # output to file
-    #------------------------------------------------------------------------------
-    print "Processor Name:",get_processor_name()
-
-    #copy script must be before screen_catcher
-
-    print 'output_dir',kwargs['output_dir']
-    
-    copy_code_files(kwargs['output_dir'],__file__, 
-             dirname(project.__file__)+sep+ project.__name__+'.py' )
-
-    store_parameters(**kwargs)
-
-    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
-
-    print "Processor Name:",get_processor_name()
-    
-    #-----------------------------------------------------------------------
-    # Domain definitions
-    #-----------------------------------------------------------------------
-
-    # Read in boundary from ordered sts file
-    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_name))
-
-    # Reading the landward defined points, this incorporates the original clipping
-    # polygon minus the 100m contour
-    landward_bounding_polygon = read_polygon(project.landward_dir)
-
-    # Combine sts polyline with landward points
-    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
-    
-    # counting segments
-    N = len(urs_bounding_polygon)-1
-
-    #boundary tags refer to project.landward 5 points equals 6 segments start at N
-    boundary_tags={'back': [N+1,N+2,N+3,N+4], 'side': [N,N+5], 'ocean': range(N)}
-
-    
-    #--------------------------------------------------------------------------
-    # 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
-    #--------------------------------------------------------------------------
-
-    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
-    # causes problems with the ability to cache set quantity which takes alot of times
-        
-    print 'start create mesh from regions'
-
-    create_mesh_from_regions(bounding_polygon,
-                         boundary_tags=boundary_tags,
-                         maximum_triangle_area=project.res_poly_all,
-                         interior_regions=project.interior_regions,
-                         filename=project.meshes_dir_name,
-                         use_cache=True,
-                         verbose=True)
-   
-    #-------------------------------------------------------------------------
-    # Setup computational domain
-    #-------------------------------------------------------------------------
-    print 'Setup computational domain'
-
-    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
-    print 'memory usage before del domain',mem_usage()
-       
-    print domain.statistics()
-    print 'triangles',len(domain)
-    
-    kwargs['act_num_trigs']=len(domain)
-
-
-    #-------------------------------------------------------------------------
-    # Setup initial conditions
-    #-------------------------------------------------------------------------
-   print 'Setup initial conditions'
-
-    #following sets the stage/water to be offcoast only
-    IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
-                             geo_reference = domain.geo_reference)
-    domain.set_quantity('stage', IC)
-    #domain.set_quantity('stage',kwargs['tide'] )
-    domain.set_quantity('friction', kwargs['friction']) 
-    
-    print 'Start Set quantity',kwargs['elevation_file']
-
-    domain.set_quantity('elevation', 
-                        filename = kwargs['elevation_file'],
-                        use_cache = False,
-                        verbose = True,
-                        alpha = kwargs['alpha'])
-    print 'Finished Set quantity'
-
-##   #------------------------------------------------------
-##    # Distribute domain to implement parallelism !!!
-##    #------------------------------------------------------
-##
-##    if numprocs > 1:
-##        domain=distribute(domain)
-
-    #------------------------------------------------------
-    # Set domain parameters
-    #------------------------------------------------------
-    print 'domain id', id(domain)
-    domain.set_name(kwargs['scenario_name'])
-    domain.set_datadir(kwargs['output_dir'])
-    domain.set_default_order(2) # Apply second order scheme
-    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
-    domain.set_store_vertices_uniquely(False)
-    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-    domain.tight_slope_limiters = 1
-    print 'domain id', id(domain)
-
-    #-------------------------------------------------------------------------
-    # Setup boundary conditions 
-    #-------------------------------------------------------------------------
-    print 'Available boundary tags', domain.get_boundary_tags()
-    print 'domain id', id(domain)
-   
-    boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
-
-    Br = Reflective_boundary(domain)
-    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-    
-    print 'Available boundary tags', domain.get_boundary_tags()
-    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
-                   domain, mean_stage= project.tide,
-                   time_thinning=1,
-                   default_boundary=Bd,
-                   use_cache=True,
-                   verbose = True,
-                   boundary_polygon=bounding_polygon)
-
-    domain.set_boundary({'back': Br,
-                         'side': Bd,
-                         'ocean': Bf}) 
-
-    kwargs['input_start_time']=domain.starttime
-
-    print'finish set boundary'
-
-    #----------------------------------------------------------------------------
-    # Evolve system through time
-    #--------------------------------------------------------------------
-    t0 = time.time()
-
-    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
-                       ,skip_initial_step = False): 
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')
-
-    x, y = domain.get_maximum_inundation_location()
-    q = domain.get_maximum_inundation_elevation()
-
-    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
-
-    print 'That took %.2f seconds' %(time.time()-t0)
-
-    #kwargs 'completed' must be added to write the final parameters to file
-    kwargs['completed']=str(time.time()-t0)
-     
-    store_parameters(**kwargs)
-     
-    print 'memory usage before del domain1',mem_usage()
-    
-    
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    
-    kwargs={}
-    kwargs['finaltime']=project.finaltime
-    kwargs['output_dir']=project.output_run_time_dir
-    kwargs['elevation_file']=project.combined_dir_name+'.pts'
-    kwargs['scenario_name']=project.scenario_name
-    kwargs['tide']=project.tide
-    kwargs['alpha'] = project.alpha
-    kwargs['friction']=project.friction
-    #kwargs['est_num_trigs']=project.trigs_min
-    #kwargs['num_cpu']=numprocs
-    #kwargs['host']=project.host
-    #kwargs['res_factor']=project.res_factor
-    #kwargs['starttime']=project.starttime
-    #kwargs['yieldstep']=project.yieldstep
-    #kwargs['file_name']=project.home+'detail.csv'
-    #kwargs['ab_time']=project.time
-    #kwargs['res_factor']= project.res_factor
-    #kwargs['user']=project.user
-    #kwargs['time_thinning'] = project.time_thinning
-    #kwargs['dir_comment']=project.dir_comment
-      
-
-    run_model(**kwargs)
-      
-   
Index: uga_work/production/Broome_2008/run_building_inundation.py
===================================================================
--- /anuga_work/production/Broome_2008/run_building_inundation.py	(revision 7159)
+++ 	(revision )
@@ -1,38 +1,0 @@
-""" Script for calculating maximum inundation at building locations
-
-    Inputs:
-
-    swwfile             - name of sww file
-                        - assume that all conserved quantities have been stored
-                    
-    buildings_filename  - name of file containing building data, sourced
-                          from NBED
-
-    Output:
-    
-    - building_filename (.csv) augmented to include maximum depth, momentum and
-      velocity and stored in same location as swwfile
-      Name = augmented_buildings.csv
-      
-    """
-
-from os import sep
-
-from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
-import project
-
-# Inputs
-#timestampdir = '20060704_063005' # HAT
-#timestampdir = '20060704_063112' # LAT
-#timestampdir = '20060704_063234' # MSL
-timestampdir = '20060515_001733' # DTED MSL
-file_loc = project.outputdir + timestampdir + sep 
-swwfile = file_loc + project.basename + '.sww'
-buildings_filename = project.buildings_filename
-buildings_filename_out = file_loc + project.buildings_filename_out
-
-#add_depth_and_momentum2csv(swwfile, buildings_filename, buildings_filename_out)
-inundation_damage(swwfile, buildings_filename, buildings_filename_out)
-
-print '\n Augmented building file written to %s \n' \
-      %buildings_filename_out
