Index: /anuga_work/production/hobart_2006/make_report.py
===================================================================
--- /anuga_work/production/hobart_2006/make_report.py	(revision 3731)
+++ /anuga_work/production/hobart_2006/make_report.py	(revision 3732)
@@ -69,5 +69,5 @@
     
 # User defined inputs
-report_title = 'Tsunami impact modelling for the Tasmania: %s' %scenario_name.title()
+report_title = 'Tsunami impact modelling for Tasmania: %s' %scenario_name.title()
 
 # WA DLI data
Index: /anuga_work/production/hobart_2006/project.py
===================================================================
--- /anuga_work/production/hobart_2006/project.py	(revision 3731)
+++ /anuga_work/production/hobart_2006/project.py	(revision 3732)
@@ -169,4 +169,10 @@
 northingmax25_2 = 5212052.309
 
+# overall 25m clipping region - covering Hobart and Bruny
+eastingmin25_3 = 523672.502
+eastingmax25_3 = 554867.24
+northingmin25_3 = 5203120.035
+northingmax25_3 = 5258511.857
+
 ###############################
 # Interior region definitions
Index: /anuga_work/production/hobart_2006/report/data.tex
===================================================================
--- /anuga_work/production/hobart_2006/report/data.tex	(revision 3731)
+++ /anuga_work/production/hobart_2006/report/data.tex	(revision 3732)
@@ -23,8 +23,10 @@
 10 and 20m contour lines.  
 
-With respect to the offshore data, the Hobart Port Authority provided state digital 
+With respect to the offshore data, the Tasmanian Ports Corporation Pty Ltd (Tasports) 
+provided state digital 
 fairsheet data around
 Hobart, in particular the Derwent River. 
-This data covers a very small geographic area. Additional data was provided from
+This data covers a very small geographic area but high density. 
+Additional data was provided from
 University of Tasmania.
 The Australian Hydrographic Office (AHO) has supplied extensive
@@ -37,5 +39,5 @@
 
 %Figure \ref{fig:contours_compare}(a) shows the contour lines for
-%HAT, MSL and LAT for teh study area using the data described above.
+%HAT, MSL and LAT for the study area using the data described above.
 
 
@@ -62,5 +64,7 @@
 \begin{tabular}{|l|l|}\hline
 Data & Specifications \\ \hline 
-DLI & Onshore, 12.5m DEM and orthophotography \\ \hline
+DPIW & Onshore, 12.5m DEM  \\ \hline
+UTAS & Offshore, sonar soundings  \\ \hline
+Tasports  & Offshore, sonar soundings  \\ \hline
 AHO & Offshore, fairsheet data for the study region \\ \hline
 \end{tabular}
Index: /anuga_work/production/hobart_2006/run_hobart_usepts.py
===================================================================
--- /anuga_work/production/hobart_2006/run_hobart_usepts.py	(revision 3732)
+++ /anuga_work/production/hobart_2006/run_hobart_usepts.py	(revision 3732)
@@ -0,0 +1,242 @@
+"""Script for running a tsunami inundation scenario for Hobart, TAS, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a tsunami wave generated by MOST.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
+"""
+#-------------------------------------------------------------------------------# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+from shutil import copy
+from os import mkdir, access, F_OK
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary, Time_boundary, File_boundary
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+
+# creates copy of code in output dir if dir doesn't exist
+if access(project.outputtimedir,F_OK) == 0 :
+    mkdir (project.outputtimedir)
+copy (project.codedirname, project.outputtimedir + project.codename)
+copy (project.codedir + 'run_hobart_usepts.py', project.outputtimedir + 'run_hobart_usepts.py')
+print'output dir', project.outputtimedir
+
+#normal screen output is stored in 
+screen_output_name = project.outputtimedir + "screen_output.txt"
+screen_error_name = project.outputtimedir + "screen_error.txt"
+
+#used to catch screen output to file
+sys.stdout = Screen_Catcher(screen_output_name)
+#sys.stderr = Screen_Catcher(screen_output_name)
+sys.stderr = Screen_Catcher(screen_error_name)
+
+print 'USER:    ', project.user
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+#-------------------------------------------------------------------------------
+
+# filenames
+#onshore_dem_name = project.onshore_dem_name
+onshore_offshore_dem_name = project.onshore_offshore_dem_name
+onshore_offshore_dem_name_25 = project.onshore_offshore_dem_name_25
+meshname = project.meshname+'.msh'
+source_dir = project.boundarydir
+
+copied_files = False
+
+# create DEM from 50m asc data
+convert_dem_from_ascii2netcdf(onshore_offshore_dem_name, use_cache=True, verbose=True)
+
+# creates pts file for combined 50m DEM and make a Geospatial data object
+dem2pts(onshore_offshore_dem_name, use_cache=True, verbose=True)
+G = Geospatial_data(file_name = project.onshore_offshore_dem_name + '.pts')
+
+# clip 50m pts based on interior regions - want 50m data OUTSIDE of these polygons
+clip_regions = [project.poly_hobart1, project.poly_hobart2, \
+                project.poly_hobart3, project.poly_hobart4]
+
+# set up initial value as Geospatial data object
+U = Geospatial_data(clip_regions[0])
+allG = G.clip_outside(U)
+for j in clip_regions[1:]:
+    U = Geospatial_data(j)
+    allG += G.clip_outside(U)
+
+print 'created outside clipping'
+# clip 50m to be inside bounding polygon
+allG = allG.clip(Geospatial_data(project.polyAll))
+print 'finished 50m clipping'
+
+# create DEM from asc data - 25m data
+convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True)
+
+# creates pts file for onshore DEM - 25 and make a Geospatial data object
+dem2pts(onshore_dem_name_offshore_25,
+        easting_min=project.eastingmin25_3,
+        easting_max=project.eastingmax25_3,
+        northing_min=project.northingmin25_3,
+        northing_max= project.northingmax25_3,
+        use_cache=True, verbose=True)
+G = Geospatial_data(file_name = project.onshore_offshore_dem_name_25 + '.pts')
+
+# clip 25m pts based on interior regions - want 25m data INSIDE these polygons
+for j in clip_regions[1:]:
+    allG += G.clip(j)
+
+print 'finished 25m clipping'
+
+allG.export_points_file(project.combined_dem_name_3 + '.pts')
+print 'exported points'
+'''
+print 'adding offshore data sets'
+G = Geospatial_data(file_name = project.offshore_dem_name_local1 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_local2 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_local3 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_local4 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho2 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho3 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho4 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho5 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho6 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho7 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho8 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho9 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho10 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho11 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho12 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho13 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho14 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho15 + '.xya')+\
+    Geospatial_data(file_name = project.offshore_dem_name_aho16 + '.xya')#+\
+    #eospatial_data(file_name = project.onshore_dem_name_25 + '.pts')
+G.export_points_file(project.combined_dem_name_3 + '.pts')
+'''
+#----------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with a tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+# use 75 for onshore components (12.5m DEM)
+hobart_res = 50000
+bruny_res = 50000
+interior_regions = [[project.poly_hobart1, hobart_res],
+                    [project.poly_hobart2, hobart_res],
+                    [project.poly_hobart3, hobart_res],
+                    [project.poly_hobart4, hobart_res]]
+
+print 'number of interior regions', len(interior_regions)
+
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project.polyAll,
+           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
+                              'e3': [3], 'e4':[4], 'e5': [5],
+			      'e6': [6], 'e7': [7], 'e8': [8],
+			      'e9': [9], 'e10': [10], 'e11': [11],
+			      'e12': [12], 'e13': [13], 'e14': [14],
+			      'e15': [15]},
+           'maximum_triangle_area': 250000,
+           'filename': meshname,           
+           'interior_regions': interior_regions},
+          verbose = True, evaluate=True)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+domain = Domain(meshname, use_cache = False, verbose = True)
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+print domain.statistics()
+
+domain.set_name(project.basename)
+domain.set_datadir(project.outputtimedir)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
+
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+tide = 0.0
+domain.set_quantity('stage', tide)
+domain.set_quantity('friction', 0.0) 
+domain.set_quantity('elevation', 
+                    filename = project.combined_dem_name + '.pts',
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.01
+                    )
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
+                    domain, verbose = True)
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+# 7 min square wave starting at 1 min, 6m high
+Bw = Time_boundary(domain = domain,
+                   f=lambda t: [(60<t<480)*10, 0, 0])
+
+# 
+domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
+		      'e5': Bd,  'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd,
+		      'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf,
+ 		      'e15': Bf} )
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 240, finaltime = 6800): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'e14')     
+
+for t in domain.evolve(yieldstep = 30, finaltime = 9000
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'e14')     
+
+for t in domain.evolve(yieldstep = 240, finaltime = 15000
+                       ,skip_initial_step = True): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'e14') 
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
+print 'finished'
