Index: anuga_work/production/FSM/build_kolonia.py
===================================================================
--- anuga_work/production/FSM/build_kolonia.py	(revision 5205)
+++ anuga_work/production/FSM/build_kolonia.py	(revision 5205)
@@ -0,0 +1,154 @@
+"""Script for running tsunami inundation scenario for Perth, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.output_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+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 import Dirichlet_boundary
+#from anuga.shallow_water import File_boundary
+#from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+#from anuga.pmesh.mesh_interface import create_mesh_from_regions
+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.combined_dir_name",project.combined_dir_name
+
+# topography directory filenames
+onshore_dir_name = project.onshore_dir_name
+onshore_dir_name1 = project.onshore_dir_name1
+#island_in_dir_name = project_urs.island_in_dir_name
+Multibeam_dir_name = project.Multibeam_dir_name
+added_data_dir_name =project.added_data_dir_name
+
+print'create Geospatial data1 objects from topographies',onshore_dir_name + '.txt'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.txt')
+
+print'create Geospatial data1 objects from topographies',onshore_dir_name1 + '.txt'
+G2 = Geospatial_data(file_name = onshore_dir_name1 + '.txt')
+
+print'create Geospatial data2 objects from multibeam', Multibeam_dir_name + '.txt'
+G_off1 = Geospatial_data(file_name = Multibeam_dir_name + '.txt')
+
+print'create Geospatial data objects from added data',added_data_dir_name + '.txt'
+G_off2 = Geospatial_data(file_name = added_data_dir_name + '.txt')
+
+print'add all geospatial objects'
+Gt1 = G1.clip(Geospatial_data(project.poly_west_reef))+G1.clip(Geospatial_data(project.poly_west))+G1.clip(Geospatial_data(project.poly_centerKolonia))+G1.clip(Geospatial_data(project.poly_east))
+Gt2= G2.clip(Geospatial_data(project.poly_topo))
+Gt3= G_off1.clip_outside(Geospatial_data(project.poly_west_reef)).clip_outside(Geospatial_data(project.poly_west)).clip_outside(Geospatial_data(project.poly_centerKolonia)).clip_outside(Geospatial_data(project.poly_east))
+                         
+G=Gt1+Gt2+Gt3+G_off2
+                         
+#print'clip combined geospatial object by bounding polygon'
+#G_clipped = G.clip(project_urs.poly_all)
+#FIXME: add a clip function to pts
+#print'shape of clipped data', G_clipped.get_data_points().shape
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+print'export',project.combined_dir_name+ '.txt'
+G.export_points_file(project.combined_dir_name+ '.txt')
+
+
+
+'''
+print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
+G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww
+
+print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
+print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
+
+#import sys; sys.exit()
+
+#if access(project.boundaries_dir,F_OK) == 0:
+#    mkdir (project.boundaries_dir)
+
+from caching import cache
+cache(urs2sww,
+      (boundaries_in_dir_name,
+       project.boundaries_dir_name1), 
+      {'verbose': True,
+       'minlat': project.south_boundary,
+       'maxlat': project.north_boundary,
+       'minlon': project.west_boundary,
+       'maxlon': project.east_boundary,
+#       'minlat': project.south,
+#       'maxlat': project.north,
+#       'minlon': project.west,
+#       'maxlon': project.east,
+       'mint': 0, 'maxt': 40000,
+#       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': project.tide,
+#       'zscale': 1,                 #Enhance tsunami
+       'fail_on_NaN': False},
+       verbose = True,
+       )
+#       dependencies = source_dir + project.boundary_basename + '.sww')
+
+'''
+
+
+
+
+
+
+
Index: anuga_work/production/FSM/changeformat.py
===================================================================
--- anuga_work/production/FSM/changeformat.py	(revision 5205)
+++ anuga_work/production/FSM/changeformat.py	(revision 5205)
@@ -0,0 +1,16 @@
+import os
+import sys
+os.chdir('E:/')
+polyline_file = open('Copy of LIDAR KOLONIA_REDUCED_5m.txt','r')
+lines=polyline_file.readlines()
+polyline_file.close()
+
+B=open('lidar_data_5m.txt','w')
+for line in lines[0:]:
+       field=line.split(' ')
+       x=float(field[0])
+       y=float(field[1])
+       z=float(field[2])
+       s='%f,%f,%f\n' %(x,y,z)
+       B.write(s)
+B.close()
Index: anuga_work/production/FSM/project.py
===================================================================
--- anuga_work/production/FSM/project.py	(revision 5205)
+++ anuga_work/production/FSM/project.py	(revision 5205)
@@ -0,0 +1,196 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
+#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
+from anuga.utilities.system_tools import get_user_name, get_host_name
+
+# file and system info
+#---------------------------------
+#codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
+user = get_user_name()
+host = get_host_name()
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+
+#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
+
+#Making assumptions about the location of scenario data
+state = 'sw_pacific'
+scenario_name = 'Kolonia_fixed_wave'
+scenario = 'FSM'
+
+tide =0.2.
+
+alpha = 0.1
+friction=0.01
+starttime=10000
+midtime=21600
+finaltime=500
+export_cellsize=50
+setup='final'
+source='test'
+
+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
+
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(scenario_name)+'_'+str(user)
+
+# onshore data 5m countour
+onshore_name = 'lidar_data_5m' # original'
+onshore_name1 = 'topo_25mcell'
+# AHO + DPI data + colin French coastline
+#coast_name = 'waterline'
+Multibeam_name = 'multibeam'
+added_data_name='sandedge'
+
+
+#final topo name
+combined_name ='Kolonia_combined_elevation'
+combined_smaller_name = 'Kolonia_combined_elevation_smaller'
+
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
+
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = anuga_dir+'topographies'+sep
+#topographies_time_dir = topographies_dir+build_time+sep
+
+# input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+onshore_in_dir_name1 = topographies_in_dir + onshore_name1
+Multibeam_in_dir_name = topographies_in_dir + Multibeam_name
+added_data_in_dir_name = topographies_in_dir + added_data_name
+
+onshore_dir_name = topographies_dir + onshore_name
+onshore_dir_name1 = topographies_dir + onshore_name1
+Multibeam_dir_name = topographies_dir + Multibeam_name
+added_data_dir_name = topographies_dir + added_data_name
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+#combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
+#combined_time_dir_final_name = topographies_time_dir + combined_final_name
+
+meshes_dir = anuga_dir+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+
+
+#boundaries_source = '1'
+
+if source =='dampier':
+    boundaries_name = 'broome_3854_17042007' #Dampier gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'dampier'+sep+'1_10000'+sep
+
+if source=='onslow':
+    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
+    
+if source=='exmouth':
+    boundaries_name = 'broome_3103_18052007' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'exmouth'+sep+'1_10000'+sep
+
+if source=='test':
+    boundaries_name = 'other' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep
+
+
+#boundaries locations
+boundaries_in_dir_name = boundaries_in_dir + scenario_name
+boundaries_dir = anuga_dir+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + scenario_name
+#boundaries_time_dir = anuga_dir+'boundaries'+sep+build_time+sep
+#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
+
+#output locations
+output_dir = anuga_dir+'outputs'+sep
+output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
+output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+#gauge_name = 'perth.csv'
+#gauges_dir = anuga_dir+'gauges'+sep
+#gauges_dir_name = gauges_dir + gauge_name
+
+#buildings_filename = gauges_dir + 'Perth_res_Project.csv'
+#buildings_filename_out = 'Perth_res_Project_modified.csv'
+
+###############################
+# Domain definitions
+###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+poly_all = read_polygon(polygons_dir+'extent.txt')
+res_poly_all = 5000*res_factor
+
+#refzone = 50 
+
+###############################
+# Interior region definitions
+###############################
+poly_topo=read_polygon(polygons_dir+'poly_topo.txt')
+poly_west_reef = read_polygon(polygons_dir+'poly_lidar1.txt')
+poly_west = read_polygon(polygons_dir+'poly_lidar4.txt')
+poly_centerKolonia = read_polygon(polygons_dir+'poly_lidar2.txt')
+poly_east = read_polygon(polygons_dir+'poly_lidar3.txt')
+
+kolonia_res = 5000
+fine_res = 100
+
+interior_regions = [[poly_west_reef, kolonia_res],
+                    [poly_west,kolonia_res],
+                    [poly_centerKolonia,fine_res],
+                    [poly_east,kolonia_res]]
+                    
+
+boundary_tags={'lagoon_west':[0],
+               'Reef_edge':[1,2,3,4,5,6,7,8,9,10,11,12,13,14],
+               'East&land':[15,16,17],
+               'Land':[18,19]}
+
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 403481
+eastingmax = 415706
+northingmin = 766423
+northingmax = 776418
+
+
+
+
+
Index: anuga_work/production/FSM/project_clip.py
===================================================================
--- anuga_work/production/FSM/project_clip.py	(revision 5205)
+++ anuga_work/production/FSM/project_clip.py	(revision 5205)
@@ -0,0 +1,197 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
+#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
+from anuga.utilities.system_tools import get_user_name, get_host_name
+
+# file and system info
+#---------------------------------
+#codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
+user = get_user_name()
+host = get_host_name()
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+
+#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
+
+#Making assumptions about the location of scenario data
+state = 'sw_pacific'
+scenario_name = 'Kolonia_fixed_wave'
+scenario = 'FSM'
+
+tide =0
+
+alpha = 0.1
+friction=0.01
+starttime=10000
+midtime=21600
+finaltime=10000
+export_cellsize=50
+setup='final'
+source='test'
+
+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
+
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(scenario_name)+'_'+str(user)
+
+# onshore data 5m countour
+onshore_name = 'lidar_data_5m' # original'
+onshore_name1 = 'topo_25m'
+# AHO + DPI data + colin French coastline
+#coast_name = 'waterline'
+Multibeam_name = 'multibeam'
+added_data_name='sandedge'
+
+
+#final topo name
+combined_name ='Kolonia_combined_elevation'
+combined_smaller_name = 'Kolonia_combined_elevation_smaller'
+combined_name1 ='Kolonia_combined_elevation1'
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
+
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = anuga_dir+'topographies'+sep
+#topographies_time_dir = topographies_dir+build_time+sep
+
+# input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+onshore_in_dir_name1 = topographies_in_dir + onshore_name1
+Multibeam_in_dir_name = topographies_in_dir + Multibeam_name
+added_data_in_dir_name = topographies_in_dir + added_data_name
+
+onshore_dir_name = topographies_dir + onshore_name
+onshore_dir_name1 = topographies_dir + onshore_name1
+Multibeam_dir_name = topographies_dir + Multibeam_name
+added_data_dir_name = topographies_dir + added_data_name
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_dir_name1 = topographies_dir + combined_name1
+#combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
+#combined_time_dir_final_name = topographies_time_dir + combined_final_name
+
+meshes_dir = anuga_dir+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+
+
+#boundaries_source = '1'
+
+if source =='dampier':
+    boundaries_name = 'broome_3854_17042007' #Dampier gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'dampier'+sep+'1_10000'+sep
+
+if source=='onslow':
+    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
+    
+if source=='exmouth':
+    boundaries_name = 'broome_3103_18052007' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'exmouth'+sep+'1_10000'+sep
+
+if source=='test':
+    boundaries_name = 'other' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep
+
+
+#boundaries locations
+boundaries_in_dir_name = boundaries_in_dir + scenario_name
+boundaries_dir = anuga_dir+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + scenario_name
+#boundaries_time_dir = anuga_dir+'boundaries'+sep+build_time+sep
+#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
+
+#output locations
+output_dir = anuga_dir+'outputs'+sep
+output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
+output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+#gauge_name = 'perth.csv'
+#gauges_dir = anuga_dir+'gauges'+sep
+#gauges_dir_name = gauges_dir + gauge_name
+
+#buildings_filename = gauges_dir + 'Perth_res_Project.csv'
+#buildings_filename_out = 'Perth_res_Project_modified.csv'
+
+###############################
+# Domain definitions
+###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+poly_all = read_polygon(polygons_dir+'extent.txt')
+res_poly_all = 5000*res_factor
+
+#refzone = 50 
+
+###############################
+# Interior region definitions
+###############################
+
+poly_west_reef = read_polygon(polygons_dir+'poly_lidar1.txt')
+poly_west = read_polygon(polygons_dir+'poly_lidar4.txt')
+poly_centerKolonia = read_polygon(polygons_dir+'poly_lidar2.txt')
+poly_east = read_polygon(polygons_dir+'poly_lidar3.txt')
+
+kolonia_res = 5000
+fine_res = 100
+
+interior_regions = [[poly_west_reef, kolonia_res],
+                    [poly_west,kolonia_res],
+                    [poly_centerKolonia,fine_res],
+                    [poly_east,kolonia_res]]
+                    
+
+boundary_tags={'lagoon_west':[0],
+               'Reef_edge':[1,2,3,4,5,6,7,8,9,10,11,12,13,14],
+               'East&land':[15,16,17],
+               'Land':[18,19]}
+
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# exporting asc grid
+eastingmin = 403481
+eastingmax = 415706
+northingmin = 766423
+northingmax = 776418
+
+
+
+
+
Index: anuga_work/production/FSM/run_kolonia.py
===================================================================
--- anuga_work/production/FSM/run_kolonia.py	(revision 5205)
+++ anuga_work/production/FSM/run_kolonia.py	(revision 5205)
@@ -0,0 +1,292 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated tsunami generated with URS code.
+
+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 Time_boundary
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary, Transmissive_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
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+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_api import distribute, numprocs, myid, barrier
+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 
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+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 kwargs
+
+    print 'output_dir',kwargs['output_dir']
+    if myid == 0:
+        copy_code_files(kwargs['output_dir'],__file__, 
+                 dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+        store_parameters(**kwargs)
+
+    barrier()
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+
+    # filenames
+#    meshes_dir_name = project.meshes_dir_name+'.msh'
+
+    # creates copy of code in output dir
+    print 'min triangles', project.trigs_min,
+    print 'Note: This is generally about 20% less than the final amount'
+
+    #--------------------------------------------------------------------------
+    # 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
+    #--------------------------------------------------------------------------
+    
+    if myid == 0:
+    
+        print 'start create mesh from regions'
+
+        create_mesh_from_regions(project.poly_all,
+                             boundary_tags=project.boundary_tags,
+                             maximum_triangle_area=project.res_poly_all,
+                             interior_regions=project.interior_regions,
+                             filename=project.meshes_dir_name+'.msh',
+                             use_cache=True,
+                             verbose=True)
+    barrier()
+
+    scenario='fixed_wave'
+    
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
+    #above don't work
+    domain = Domain(project.meshes_dir_name+'.msh', use_cache=True, 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
+    #-------------------------------------------------------------------------
+    if myid == 0:
+
+        print 'Setup initial conditions'
+
+        from polygon import Polygon_function
+        #following sets the stage/water to be offcoast only
+#        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
+#                                 geo_reference = domain.geo_reference)
+#        domain.set_quantity('stage', IC)
+        domain.set_quantity('stage',kwargs['tide'] )
+#        domain.set_quantity('stage', kwargs['tide'])
+        domain.set_quantity('friction', kwargs['friction']) 
+        
+        print 'Start Set quantity',kwargs['bathy_file']
+
+        domain.set_quantity('elevation', 
+                            filename = kwargs['bathy_file'],
+                            use_cache = False,
+                            verbose = True,
+                            alpha = kwargs['alpha'])
+        print 'Finished Set quantity'
+    barrier()
+    
+    #------------------------------------------------------
+    # Distribute domain to implement parallelism !!!
+    #------------------------------------------------------
+
+    if numprocs > 1:
+        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(kwargs['aa_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
+    #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+
+    #------------------------------------------------------------------------------
+    # Setup information for slide scenario (to be applied 1 min into simulation
+    #------------------------------------------------------------------------------
+
+    if scenario == 'slide':
+        # Function for submarine slide
+        from anuga.shallow_water.smf import slide_tsunami  
+        tsunami_source = slide_tsunami(length=35000.0,
+                                       depth=project.slide_depth,
+                                       slope=10.0,
+                                       thickness=800.0, 
+                                       x0=project.slide_origin[0], 
+                                       y0=project.slide_origin[1], 
+                                       alpha=0.0, 
+                                       domain=domain,
+                                       verbose=True)
+    
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
+
+    Bt = Transmissive_boundary(domain)
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([0,0,0])
+
+    if scenario == 'fixed_wave':
+    # FIXME (Ole): Change this to Transmissive Momentum as
+    # suggested by Rajaraman (ticket:208)
+        from math import sin, pi
+        
+        Bw = Time_boundary(domain = domain,
+                           f=lambda t: [(120<t<220)*3.0*sin(t*pi*2/10), 0, 0])
+        domain.set_boundary({'lagoon_west':Bt,
+                             'Land':Bt,
+                             'Reef_edge':Bw,
+                             'East&land':Bt})
+
+    # boundary conditions for slide scenario
+    if scenario == 'slide':
+        domain.set_boundary({'ocean':Bd,'land':Bd})
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(kwargs['yieldstep'], finaltime = kwargs['finaltime']): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'Reef_edge')     
+
+        
+    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)
+    
+    if myid==0:
+        store_parameters(**kwargs)
+    barrier
+    
+    print 'memory usage before del domain1',mem_usage()
+
+def export_model(**kwargs):    
+    #store_parameters(**kwargs)
+    
+#    print 'memory usage before del domain',mem_usage()
+    #del domain
+    print 'memory usage after del domain',mem_usage()
+    
+    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
+    print'swwfile',swwfile
+    
+    export_grid(swwfile, extra_name_out = 'FSM',
+            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            timestep = None,
+            reduction = max,
+            cellsize = 25,
+            NODATA_value = -1E-030,
+            easting_min = project.eastingmin,
+            easting_max = project.eastingmax,
+            northing_min = project.northingmin,
+            northing_max = project.northingmax,
+            verbose = False,
+            origin = None,
+            datum = 'WGS84',
+            format = 'txt')
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    kwargs={}
+    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['midtime']=project.midtime
+    kwargs['finaltime']=project.finaltime
+    kwargs['output_dir']=project.output_run_time_dir
+    kwargs['bathy_file']=project.combined_dir_name+'.txt'
+#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
+    kwargs['file_name']=project.home+'detail.csv'
+    kwargs['aa_scenario_name']=project.scenario_name
+    kwargs['ab_time']=project.time
+    kwargs['res_factor']= project.res_factor
+    kwargs['tide']=project.tide
+    kwargs['user']=project.user
+    kwargs['alpha'] = project.alpha
+    kwargs['friction']=project.friction
+    kwargs['time_thinning'] = project.time_thinning
+    kwargs['dir_comment']=project.dir_comment
+    kwargs['export_cellsize']=project.export_cellsize
+    
+
+    run_model(**kwargs)
+      
+    if myid==0:
+        export_model(**kwargs)
+    barrier
Index: anuga_work/production/FSM/test_clip
===================================================================
--- anuga_work/production/FSM/test_clip	(revision 5205)
+++ anuga_work/production/FSM/test_clip	(revision 5205)
@@ -0,0 +1,103 @@
+# 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.combined_dir_name",project.combined_dir_name
+
+# topography directory filenames
+onshore_dir_name = project.onshore_dir_name
+onshore_dir_name1 = project.onshore_dir_name
+#island_in_dir_name = project_urs.island_in_dir_name
+Multibeam_dir_name = project.Multibeam_dir_name
+added_data_dir_name =project.added_data_dir_name
+
+print'create Geospatial data1 objects from topographies',onshore_dir_name + '.txt'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.txt')
+
+print'create Geospatial data1 objects from topographies',onshore_dir_name1 + '.txt'
+G2 = Geospatial_data(file_name = onshore_dir_name1 + '.txt')
+
+print'create Geospatial data2 objects from multibeam', Multibeam_dir_name + '.txt'
+G_off1 = Geospatial_data(file_name = Multibeam_dir_name + '.txt')
+
+print'create Geospatial data objects from added data',added_data_dir_name + '.txt'
+G_off2 = Geospatial_data(file_name = added_data_dir_name + '.txt')
+
+print'add all geospatial objects'
+Gt1 = G1.clip(Geospatial_data(project.poly_west_reef))+G1.clip(Geospatial_data(project.poly_west))+G1.clip(Geospatial_data(project.poly_centerKolonia))+G1.clip(Geospatial_data(project.poly_east))
+Gt2= G2.clip_outside(Geospatial_data(project.poly_centerKolonia)).clip_outside(Geospatial_data(project.poly_east))
+Gt3= G_off1.clip_outside(Geospatial_data(project.poly_west_reef)).clip_outside(Geospatial_data(project.poly_west)).clip_outside(Geospatial_data(project.poly_centerKolonia)).clip_outside(Geospatial_data(project.poly_east))
+                         
+G=Gt1+Gt2+Gt3+G_off2
+                         
+#print'clip combined geospatial object by bounding polygon'
+#G_clipped = G.clip(project_urs.poly_all)
+#FIXME: add a clip function to pts
+#print'shape of clipped data', G_clipped.get_data_points().shape
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+print'export',project.combined_dir_name+ '.txt'
+G.export_points_file(project.combined_dir_name+ '.txt')
+
+
+
+'''
+print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
+G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww
+
+print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
+print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
+
+#import sys; sys.exit()
+
+#if access(project.boundaries_dir,F_OK) == 0:
+#    mkdir (project.boundaries_dir)
+
+from caching import cache
+cache(urs2sww,
+      (boundaries_in_dir_name,
+       project.boundaries_dir_name1), 
+      {'verbose': True,
+       'minlat': project.south_boundary,
+       'maxlat': project.north_boundary,
+       'minlon': project.west_boundary,
+       'maxlon': project.east_boundary,
+#       'minlat': project.south,
+#       'maxlat': project.north,
+#       'minlon': project.west,
+#       'maxlon': project.east,
+       'mint': 0, 'maxt': 40000,
+#       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': project.tide,
+#       'zscale': 1,                 #Enhance tsunami
+       'fail_on_NaN': False},
+       verbose = True,
+       )
+#       dependencies = source_dir + project.boundary_basename + '.sww')
+
+'''
+
+
+
+
+
+
+
+
Index: anuga_work/production/FSM/test_clip.py
===================================================================
--- anuga_work/production/FSM/test_clip.py	(revision 5205)
+++ anuga_work/production/FSM/test_clip.py	(revision 5205)
@@ -0,0 +1,137 @@
+"""Script for running tsunami inundation scenario for Perth, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_clip.py
+The output sww file is stored in project_clip.output_time_dir 
+
+The scenario is defined by a triangular mesh created from project_clip.polygon,
+the elevation data and a simulated submarine landslide.
+
+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 import Dirichlet_boundary
+#from anuga.shallow_water import File_boundary
+#from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+#from anuga.pmesh.mesh_interface import create_mesh_from_regions
+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_clip   # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+copy_code_files(project_clip.output_build_time_dir,__file__, 
+               dirname(project_clip.__file__)+sep+ project_clip.__name__+'.py' )
+
+start_screen_catcher(project_clip.output_build_time_dir)
+
+print 'USER:    ', project_clip.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_clip.combined_dir_name",project_clip.combined_dir_name
+
+# topography directory filenames
+total=project_clip.combined_dir_name
+
+print'create Geospatial data1 objects from topographies',total + '.txt'
+G1 = Geospatial_data(file_name = total + '.txt')
+
+print'add all geospatial objects'
+
+Gt= G1.clip_outside(Geospatial_data(project_clip.poly_west_reef)).clip_outside(Geospatial_data(project_clip.poly_west)).clip_outside(Geospatial_data(project_clip.poly_centerKolonia)).clip_outside(Geospatial_data(project_clip.poly_east))
+                         
+#print'clip combined geospatial object by bounding polygon'
+#G_clipped = G.clip(project_clip_urs.poly_all)
+#FIXME: add a clip function to pts
+#print'shape of clipped data', G_clipped.get_data_points().shape
+
+print'export combined DEM file'
+if access(project_clip.topographies_dir,F_OK) == 0:
+    mkdir (project_clip.topographies_dir)
+print'export',project_clip.combined_dir_name+ '.txt'
+Gt.export_points_file(project_clip.combined_dir_name1+ '.txt')
+
+
+
+'''
+print'project_clip.combined_dir_name + 1.xya',project_clip.combined_dir_name + '1.xya'
+G_all=Geospatial_data(file_name = project_clip.combined_dir_name + '1.xya')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project_clip.combined_dir_name+'_small1' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project_clip.combined_dir_name+'_other1' + '.xya')
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project_clip.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww
+
+print 'minlat=project_clip.north_boundary, maxlat=project_clip.south_boundary',project_clip.north_boundary, project_clip.south_boundary
+print 'minlon= project_clip.west_boundary, maxlon=project_clip.east_boundary',project_clip.west_boundary, project_clip.east_boundary
+
+#import sys; sys.exit()
+
+#if access(project_clip.boundaries_dir,F_OK) == 0:
+#    mkdir (project_clip.boundaries_dir)
+
+from caching import cache
+cache(urs2sww,
+      (boundaries_in_dir_name,
+       project_clip.boundaries_dir_name1), 
+      {'verbose': True,
+       'minlat': project_clip.south_boundary,
+       'maxlat': project_clip.north_boundary,
+       'minlon': project_clip.west_boundary,
+       'maxlon': project_clip.east_boundary,
+#       'minlat': project_clip.south,
+#       'maxlat': project_clip.north,
+#       'minlon': project_clip.west,
+#       'maxlon': project_clip.east,
+       'mint': 0, 'maxt': 40000,
+#       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': project_clip.tide,
+#       'zscale': 1,                 #Enhance tsunami
+       'fail_on_NaN': False},
+       verbose = True,
+       )
+#       dependencies = source_dir + project_clip.boundary_basename + '.sww')
+
+'''
+
+
+
+
+
+
+
+
