Index: /anuga_work/development/boxingday_2007/build_boxing.py
===================================================================
--- /anuga_work/development/boxingday_2007/build_boxing.py	(revision 4297)
+++ /anuga_work/development/boxingday_2007/build_boxing.py	(revision 4297)
@@ -0,0 +1,93 @@
+"""Script for running tsunami inundation scenario for boxing day Phuket,
+Thailand. 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.
+
+John Jakeman - ANU 2007
+"""
+
+#use _read_asc() in data_manager
+
+
+#------------------------------------------------------------------------------
+# 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
+from time import localtime, strftime
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts, ferret2sww
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+from anuga.utilities.polygon import polygon_area
+
+# Application specific imports
+import project   # Definition of file names and polygons
+
+def convert_arcgis_latlon_list_to_utm(points):
+     #Used because arc gis produced csv files put lat lon in
+     #reverse order to those accpeted by convert_latlon_to_utm()
+     
+     from anuga.coordinate_transforms.geo_reference import Geo_reference
+     from anuga.coordinate_transforms.redfearn import redfearn
+     old_geo = Geo_reference() 
+     utm_points = []
+     for point in points:
+         zone, easting, northing = redfearn(float(point[1]),float(point[0]))
+         new_geo = Geo_reference(zone)
+         old_geo.reconcile_zones(new_geo) 
+         utm_points.append([easting,northing])
+
+     return utm_points, old_geo.get_zone()
+"""
+#------------------------------------------------------------------------------
+# 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.bounding_polygon_latlon",project.bounding_polygon_latlon
+
+# topography and bathymetry directory filenames
+bathymetry_file = project.bathymetry
+
+#gets time for new dir
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+# creates DEM from asc data
+print "creates DEMs from asc data"
+convert_dem_from_ascii2netcdf(bathymetry_file, basename_out=project.bathymetry_points, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(project.bathymetry_points,
+        use_cache=True, 
+        verbose=True)
+
+print'create Geospatial data objects from topographies'
+G = Geospatial_data(file_name = project.bathymetry_points + '.pts')
+
+bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon)
+print bounding_polygon
+print 'zone of bounding polygon', zone
+print 'Area of bounding polygon', polygon_area(bounding_polygon)/1000000.0
+
+print'clip combined geospatial object by bounding polygon'
+G_clipped = G.clip(bounding_polygon)
+
+print'export combined DEM file'
+G_clipped.export_points_file(project.combined_dir_name + '.pts')
+"""
+
Index: /anuga_work/development/boxingday_2007/convert.py
===================================================================
--- /anuga_work/development/boxingday_2007/convert.py	(revision 4297)
+++ /anuga_work/development/boxingday_2007/convert.py	(revision 4297)
@@ -0,0 +1,30 @@
+from anuga.geospatial_data.geospatial_data import *
+from anuga.coordinate_transforms.redfearn import *
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts,xya2pts
+
+#xya2pts('thaicoas_9', verbose=True)
+print'create G1'
+G1 = Geospatial_data(file_name = 'thaicoas_9' + '.xya')
+G = G1
+G.export_points_file('boxingday' + '.pts')
+
+
+# Define extent and print into csv file
+points = [(6.0,97.5),(10.0,97.5),(10.0,100.5),(6.0,100.5)]
+points, zone = convert_from_latlon_to_utm(points)
+#print 'points', points
+#print 'zone', zone
+
+file_name = 'extent' +'.csv'
+f = file(file_name, 'w')
+for i in range(len(points)):
+    #print '%f, %f' %(points[i][0],points[i][1])
+    f.write('%f, %f\n' %(points[i][0],points[i][1]))
+f.close
+
+# Define Interior regions and print into corresponding csv files
+# Coast
+
+# Island 1
+
+# Island 2
Index: /anuga_work/development/boxingday_2007/north_island.csv
===================================================================
--- /anuga_work/development/boxingday_2007/north_island.csv	(revision 4297)
+++ /anuga_work/development/boxingday_2007/north_island.csv	(revision 4297)
@@ -0,0 +1,4 @@
+97.61726334710,8.68855366473
+97.67860649750,8.68855366473
+97.67685383610,8.44142840154
+97.61200536270,8.44318106298
Index: /anuga_work/development/boxingday_2007/project.py
===================================================================
--- /anuga_work/development/boxingday_2007/project.py	(revision 4297)
+++ /anuga_work/development/boxingday_2007/project.py	(revision 4297)
@@ -0,0 +1,52 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, \
+                                    polygon_area, is_inside_polygon
+from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+from time import localtime, strftime
+
+###############################
+# Domain definitions
+###############################
+
+# Bathymetry and topography filenames
+bathymetry = 'phuket_utm47'
+
+# bounding polygon for study area
+bounding_polygon_latlon = read_polygon('phuket_boundary.csv')
+# only works for zone 47
+south = 1000
+east = -1000
+west = 1000
+north = -1000
+for points in bounding_polygon_latlon:
+    south = min(points[1],south)
+    north = max(points[1],north)
+    east = max(points[0],east)
+    west = min(points[0],west)
+
+print "south, north, east, west", south, north, east, west
+
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+#combined_dir_name = 'phuket_bathymetry'+time
+combined_dir_name = 'phuket_bathymetry'
+meshname = 'phuket_mesh'+time
+mesh_elevname = 'phuket_mesh'+time+'_elev'
+bathymetry_points = 'bathymetry_points'+time
+boundary_most_in = 'out'
+boundary_most_out = 'most_boundary_condition'
+###############################
+# Interior region definitions
+###############################
+
+# interior polygons
+poly_island1 = read_polygon('island1.csv')
+poly_north_island = read_polygon('north_island.csv')
+poly_50m_contour = read_polygon('phuket_50m_contour_small.csv')
+
+gauge_filename = 'gauges.csv'
Index: /anuga_work/development/boxingday_2007/run_boxingday.py
===================================================================
--- /anuga_work/development/boxingday_2007/run_boxingday.py	(revision 4297)
+++ /anuga_work/development/boxingday_2007/run_boxingday.py	(revision 4297)
@@ -0,0 +1,287 @@
+"""Script for running a tsunami inundation scenario for Cairns, QLD 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 directory named after the scenario, i.e
+slide or fixed_wave.
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a tsunami wave generated by a submarine mass failure.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and
+Nick Bartzis, GA - 2006
+"""
+
+def convert_arcgis_latlon_list_to_utm(points):
+     #Used because arc gis produced csv files put lat lon in
+     #reverse order to those accpeted by convert_latlon_to_utm()
+     
+     from anuga.coordinate_transforms.geo_reference import Geo_reference
+     from anuga.coordinate_transforms.redfearn import redfearn
+     old_geo = Geo_reference() 
+     utm_points = []
+     for point in points:
+         zone, easting, northing = redfearn(float(point[1]),float(point[0]))
+         new_geo = Geo_reference(zone)
+         old_geo.reconcile_zones(new_geo) 
+         utm_points.append([easting,northing])
+
+     return utm_points, old_geo.get_zone()
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import Time_boundary
+from anuga.shallow_water import File_boundary
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww
+from anuga.shallow_water.data_manager import dem2pts
+from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+#from anuga.fit_interpolate.fit import fit_to_mesh_file (does not exist)
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+
+#------------------------------------------------------------------------------
+# Define scenario as either slide or fixed_wave.
+#------------------------------------------------------------------------------
+
+#scenario = 'coseismic'
+scenario = 'fixed_wave'
+
+if os.access(scenario, os.F_OK) == 0:
+    os.mkdir(scenario)
+basename = scenario + 'source'
+
+
+#------------------------------------------------------------------------------
+# Preparation of topographic data
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+#------------------------------------------------------------------------------
+
+# Filenames
+dem_name = 'boxingday' 
+
+# Create DEM from asc data
+#convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
+# Create pts file for onshore DEM
+#dem2pts(dem_name, use_cache=True, verbose=True)
+
+
+#------------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with a tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#------------------------------------------------------------------------------
+
+remainder_res = 10000000
+islands_res = 100000
+cairns_res = 100000
+shallow_res = 500000
+#interior_regions = [[project.poly_cairns, cairns_res]#,
+#                    [project.poly_island0, islands_res],
+#                    [project.poly_island1, islands_res],
+#                    [project.poly_island2, islands_res],
+#                    [project.poly_island3, islands_res],
+#                    [project.poly_shallow, shallow_res]]
+
+# filenames
+meshname = project.meshname + '.tsh'
+mesh_elevname = project.mesh_elevname  + '.tsh'
+
+bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon)
+
+from caching import cache
+print 'start create mesh from regions'
+
+if scenario == 'coseismic':
+     _ = cache(create_mesh_from_regions,
+               bounding_polygon,
+               {'boundary_tags': {'one': [0], 'two': [1],
+                                  'three': [2], 'four': [3],
+                                  'five': [4], 'six': [5],
+                                  'seven': [6], 'eight': [7],
+                                  'nine': [8], 'ten': [9],
+                                  'eleven': [10], 'twelve': [11],
+                                  'thirteen': [12], 'fourteen': [13]},
+                'maximum_triangle_area': 5000000,
+                'filename': meshname},#,           
+               #'interior_regions': interior_regions},
+               verbose = True
+               #, evaluate=True
+               )
+
+if scenario == 'fixed_wave':
+     south = 7.40911081272
+     north = 8.71484358635
+     east = 99.1683687224
+     west = 97.513856322
+     points = [[south,west],[south,east],[north,east],[north,west]]
+
+     bounding_polygon,zone=convert_from_latlon_to_utm(points)
+
+     _ = cache(create_mesh_from_regions,
+               bounding_polygon,
+               {'boundary_tags': {'ocean_west': [0], 'bottom': [1],
+               'onshore': [2], 'top': [3]},
+                'maximum_triangle_area': 500000000,
+                'filename': meshname},#,           
+               #'interior_regions': interior_regions},
+               verbose = True
+               #, evaluate=True
+               )
+
+## cache(fit_to_mesh_file,(meshname,
+##       project.combined_dir_name + '.pts',
+##       mesh_elevname),
+##       #{'verbose': True},
+##       verbose = False
+##      )
+
+#------------------------------------------------------------------------------
+# 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(basename) 
+domain.set_datadir(scenario)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def elevation(x,y):
+     return -10.0
+     
+tide = 0.0
+domain.set_quantity('stage', tide)
+domain.set_quantity('friction', 0.0) 
+
+if scenario == 'fixed_wave':
+     # test with coarser bathymetry                        
+     domain.set_quantity('elevation',
+                         #filename = 'thaicoas_9.pts',
+                         elevation,
+                         use_cache=True,
+                         verbose=True,
+                         alpha=0.1)
+
+if scenario == 'coseismic':
+     domain.set_quantity('elevation', 
+                         filename=project.combined_dir_name + '.pts',
+                         use_cache=True,
+                         verbose=True,
+                         alpha=0.1)
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+print 'start ferret2sww'
+south = project.south 
+north = project.north
+west = project.west
+east = project.east
+
+"""
+tide = 0.0
+# Note only need to do when an SWW file for the MOST boundary doesn't exist
+cache(ferret2sww,
+      (project.boundary_most_in,
+       project.boundary_most_out), 
+      {'verbose': True,
+       'minlat': south,
+       'maxlat': north,
+       'minlon': west,
+       'maxlon': east,
+#       'origin': project.mesh_origin,
+       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': tide,
+       'zscale': 10,                 #Enhance tsunami (Does this affect result)
+       'fail_on_NaN': False,
+       'inverted_bathymetry': True},
+       #evaluate = True,
+       verbose = True)#,
+      #dependencies =  'most_results' + '.sww')
+
+"""
+print 'Available boundary tags', domain.get_boundary_tags()
+
+#Bf = File_boundary('most_results'+'.sww', 
+#                    domain, verbose = True)
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+Bw = Time_boundary(domain = domain,
+                       f=lambda t: [(60<t<3660)*50, 0, 0])
+    
+# 60 min square wave starting at 1 min, 50m high
+if scenario == 'fixed_wave':
+    domain.set_boundary({'ocean_west': Br,
+                         'bottom': Br,
+                         'onshore': Br,
+                         'top': Br})
+
+# boundary conditions for slide scenario
+if scenario == 'coseismic':
+    domain.set_boundary( {'one': Bf,'two': Bf, 'three':  Bf, 'four': Bf,
+                      'five': Bf, 'six': Bd, 'seven': Bd, 'eight': Bf,
+                      'nine': Bf, 'ten':  Bd, 'eleven': Bd, 'twelve': Bf,
+                      'thireteen': Bf, 'fourteen': Bf})
+    
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+#domain.visualise = True
+#domain.visualise_color_stage = True
+
+import time
+t0 = time.time()
+
+from Numeric import allclose
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+if scenario == 'cosesimic':
+    
+     for t in domain.evolve(yieldstep = 10, finaltime = 60):
+          print domain.quantities['elevation'].vertex_values
+          domain.write_time() 
+
+     for t in domain.evolve(yieldstep = 600, finaltime = 21600, 
+                            skip_initial_step = True):
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean_west')
+
+if scenario == 'fixed_wave':
+
+    # save every two mins leading up to wave approaching land
+    for t in domain.evolve(yieldstep = 10, finaltime = 50):
+         print domain.quantities['elevation'].vertex_values
+         print domain.quantities['stage'].vertex_values
+         
+         domain.write_time()
+
+    #save every 30 secs as wave starts inundating ashore
+    #for t in domain.evolve(yieldstep = 10, finaltime = 10000, 
+    #                       skip_initial_step = True):
+    #    domain.write_time()
+            
+print 'That took %.2f seconds' %(time.time()-t0)
