Index: /anuga_work/development/PNG_slide/export_results.py
===================================================================
--- /anuga_work/development/PNG_slide/export_results.py	(revision 5634)
+++ /anuga_work/development/PNG_slide/export_results.py	(revision 5634)
@@ -0,0 +1,76 @@
+import project_slide, os
+import sys
+
+from anuga.shallow_water.data_manager import sww2dem
+from os import sep
+
+#time_dir = '20061211_071516' #bulli
+#time_dir = '20061212_012705' #shovel
+#time_dir = '20061212_064807' #yacaaba
+# correct depths
+#time_dir = '20070121_213509' #shovel
+#time_dir = '20070119_055222' #yacaaba
+#time_dir = '20070121_213444' #bulli
+#corrected yacaaba and birubi
+#time_dir = '20070322_035214' # yacaaba
+time_dir = '20070322_035135' # birubi
+
+directory = project_slide.outputdir
+name = directory + time_dir + sep + project_slide.basename
+#name = 'test500000.035000'
+#name = 'test'
+
+is_parallel = False
+if is_parallel == True: nodes = 4
+print 'output dir:', name
+
+which_var = 2
+if which_var == 0:  # Stage
+    outname = name + '_stage-eq-thurs-test'
+    quantityname = 'stage'
+
+if which_var == 1:  # Absolute Momentum
+    outname = name + '_momentum_i1'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  
+
+if which_var == 2:  # Depth
+    outname = name + '_depth'
+    quantityname = 'stage-elevation'  
+
+if which_var == 3:  # Speed
+    outname = name + '_speed_i0'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'  #Speed
+
+if which_var == 4:  # Elevation
+    outname = name + '_elevation'
+    quantityname = 'elevation'  #Elevation
+
+if is_parallel == True:
+    for i in range(0,nodes):
+        namei = name + '_P%d_%d' %(i,nodes)
+        outnamei = outname + '_P%d_%d' %(i,nodes)
+        print 'start sww2dem for sww file %d' %(i)
+        sww2dem(namei, basename_out = outnamei,
+                    quantity = quantityname,
+                    timestep = 1,
+                    cellsize = 100,      
+                    #easting_min = project.e_min_area,
+                    #easting_max = project.e_max_area,
+                    #northing_min = project.n_min_area,
+                    #northing_max = project.n_max_area,        
+                    reduction = max, 
+                    verbose = True,
+                    format = 'asc')
+else:
+    print 'start sww2dem'
+    sww2dem(name, basename_out = outname,
+                quantity = quantityname,
+                timestep = 0,
+                cellsize = 50,      
+                easting_min = project_slide.eastingmin,
+                easting_max = project_slide.eastingmax,
+                northing_min = project_slide.northingmin,
+                northing_max = project_slide.northingmax,        
+                reduction = max, 
+                verbose = True,
+                format = 'asc')
Index: /anuga_work/development/PNG_slide/project_slide.py
===================================================================
--- /anuga_work/development/PNG_slide/project_slide.py	(revision 5634)
+++ /anuga_work/development/PNG_slide/project_slide.py	(revision 5634)
@@ -0,0 +1,209 @@
+# -*- 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, polygon_area, is_inside_polygon
+
+if sys.platform == 'win32':
+    home = getenv('INUNDATIONHOME')
+    user = getenv('USERPROFILE')
+
+else:    
+    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
+    user = getenv('LOGNAME')
+    print 'USER:', user
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
+
+#Making assumptions about the location of scenario data
+state = 'new_south_wales'
+scenario_dir_name = 'sydney_tsunami_scenario_2006'
+
+# onshore data provided by NSW LPI, Krishna merged with existing
+# PMD data to make 100m DEM
+on_offshore_name25 = 'bathy_dem25'
+on_offshore_name100 = 'bathyland100'
+nsw100 = 'dem_utm'
+
+# survey data
+offshore_name1 = 'surveyAreaA'
+offshore_name2 = 'surveyAreaB'
+offshore_name3 = 'surveyAreaC'
+
+# AHO data
+offshore_name4 = '1000003611export'
+offshore_name5 = '1000003613export'
+offshore_name6 = '1000003614export'
+offshore_name7 = '1000003627export'
+offshore_name8 = '1000003628export'
+offshore_name9 = 'AHDexport'
+
+#swollen/ all data output
+basename = 'source'
+codename = 'project_slide.py'
+
+#Derive subdirectories and filenames
+local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
+meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
+datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
+gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
+outputtimedir = outputdir + local_time + sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+
+gauge_filename = gaugedir + 'sydney_slide_gauges.csv'
+codedir = getcwd()+sep                            
+codedirname = codedir + 'project_slide.py'
+meshname = outputtimedir + 'mesh_' + basename
+
+# Necessary if using point datasets, rather than grid
+on_offshore25_dem_name = datadir + on_offshore_name25
+on_offshore100_dem_name = datadir + on_offshore_name100
+nsw_dem_name = datadir + nsw100
+offshore_dem_name1 = datadir + offshore_name1
+offshore_dem_name2 = datadir + offshore_name2
+offshore_dem_name3 = datadir + offshore_name3
+offshore_dem_name4 = datadir + offshore_name4
+offshore_dem_name5 = datadir + offshore_name5
+offshore_dem_name6 = datadir + offshore_name6
+offshore_dem_name7 = datadir + offshore_name7
+offshore_dem_name8 = datadir + offshore_name8
+offshore_dem_name9 = datadir + offshore_name9
+combined_dem_name   = datadir + 'sydney_combined_elevation'
+
+###############################
+# Domain definitions
+###############################
+
+# bounding polygon for study area
+polyAll = read_polygon(polygondir+'extent3.csv')
+
+print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
+
+###############################
+# Interior region definitions
+###############################
+
+# sydney digitized polygons
+poly_coast = read_polygon(polygondir+'coast2.csv')
+poly_syd1 = read_polygon(polygondir+'local.csv')
+
+print 'Area of local polygon', polygon_area(poly_syd1)/1000000.0
+print 'Area of coastal polygon', polygon_area(poly_coast)/1000000.0
+
+###################################################################
+# Clipping regions for export to asc and regions for clipping data
+###################################################################
+
+# clipping for 25m grid for conversion
+eastingmin25 = 324825
+eastingmax25 = 347650
+northingmin25 = 6231960
+northingmax25 = 6272543
+
+poly_25mclip = read_polygon(polygondir+'25mclip.csv')
+
+# clipping for 100m grid for conversion
+eastingmin100 = 315870
+eastingmax100 = 361607
+northingmin100 = 6227588
+northingmax100 = 6285975
+
+poly_surveyclip = read_polygon(polygondir+'surveyclip.csv')
+
+# clipping for nsw grid for conversion
+eastingmin_nsw = 315870
+eastingmax_nsw = 390480
+northingmin_nsw = 6161100 
+northingmax_nsw = 6330000
+
+poly_origsyd = [[eastingmin100,northingmin100],[eastingmax100,northingmin100],\
+                [eastingmax100,northingmax100],[eastingmin100,northingmax100]]
+
+# exporting asc grid
+eastingmin = 333000
+eastingmax = 354000
+northingmin = 6244800
+northingmax = 6274550
+
+# exporting asc grid
+eastingmin = 318922+10000
+eastingmax = 318922+45000
+northingmin = 6231749+7500
+northingmax = 6231749+48000
+
+#eastingmin = -40
+#eastingmax = 40
+#northingmin = -20
+#northingmax = 20
+
+###################################################################
+# Slide characteristics
+###################################################################
+
+# historical slides
+slide_origin_bulli = [370110,6189489]
+slide_origin_shovel = [354753,6187397]
+#slide_origin_yacaaba = [463012,6334735] original work was actually Birubi
+slide_origin_yacaaba = [476465, 6362424]
+slide_origin_birubi = [462267,6335100]
+
+# potential slides: these correspond approximately to modelled depths
+slide_origin_bulli_b = [405617,6249852]
+slide_origin_bulli_c0 = [438583,6308510]
+slide_origin_bulli_c1 = [465283,6335210]
+slide_origin_bulli_c2 = [483783,6358810]
+slide_origin_shovel_b = [396262,6252144]
+slide_origin_shovel_c0 = [431683,6311110]
+slide_origin_shovel_c1 = [458694,6340352]
+slide_origin_shovel_c2 = [475589,6363030]
+slide_origin_yacaaba_a = [359796,6190943]
+slide_origin_yacaaba_b = [397330,6251857]
+slide_origin_yacaaba_c0 = [431283,6312110]
+slide_origin_yacaaba_c1 = [459983,6339110]
+slide_origin_birubi_a = [362811,6191183]
+slide_origin_birubi_b = [403573,6250268]
+slide_origin_birubi_c0 = [437183,6309410]
+slide_origin_birubi_c2 = [483392,6359044]
+
+#bulli_depth = 2087.0
+bulli_depth = 1470.0
+bulli_length = 16840.0
+bulli_thickness = 424.0
+bulli_width = 8860.0
+bulli_density = 1.46
+bulli_slope = 4.0
+bulli_alpha = 126.0 - 90.0
+
+#shovel_depth = 968.0
+shovel_depth = 877.0
+shovel_length = 13500.0
+shovel_thickness = 165.0
+shovel_width = 4350.0
+shovel_density = 1.49
+shovel_slope = 4.0
+shovel_alpha = 118.0 - 90.0
+
+#yacaaba_depth = 1119.0
+yacaaba_depth = 938.
+yacaaba_length = 4189.
+yacaaba_thickness = 53.
+yacaaba_width = 2898.
+yacaaba_density = 1.48
+yacaaba_slope = 2.3
+yacaaba_alpha = 133.0 - 90.0
+ 
+birubi_depth = 1320.
+birubi_length = 9903.0
+birubi_width = 4150.
+birubi_density = 1.48
+birubi_slope = 3.7
+birubi_alpha = 133.0 - 90.0
+birubi_thickness = 140.
Index: /anuga_work/development/PNG_slide/run_png_test.py
===================================================================
--- /anuga_work/development/PNG_slide/run_png_test.py	(revision 5634)
+++ /anuga_work/development/PNG_slide/run_png_test.py	(revision 5634)
@@ -0,0 +1,208 @@
+"""Script for running a tsunami inundation scenario for Aitape, PNG.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_slide.py
+The output sww file is stored in project_slide.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project_slide.polygon,
+the elevation data and a tsunami wave generated by s submarine mass failure.
+
+Bridgette Lewis, GA - 2008
+"""
+
+#-------------------------------------------------------------------------------
+# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+from shutil import copy
+from os.path import dirname, basename
+from os import mkdir, access, F_OK, sep
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+
+# Application specific imports
+import project_slide              # 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 
+copy_code_files(project_slide.outputtimedir,__file__,dirname(project_slide.__file__)+sep+ project_slide.__name__+'.py' )
+myid = 0
+numprocs = 1
+start_screen_catcher(project_slide.outputtimedir, myid, numprocs)
+
+print 'USER:    ', project_slide.user
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+#-------------------------------------------------------------------------------
+
+# filenames
+on_offshore25_dem_name = project_slide.on_offshore25_dem_name
+on_offshore100_dem_name = project_slide.on_offshore100_dem_name
+nsw_dem_name = project_slide.nsw_dem_name
+meshname = project_slide.meshname+'.msh'
+
+# creates DEM from asc data
+convert_dem_from_ascii2netcdf(on_offshore25_dem_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(on_offshore100_dem_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(nsw_dem_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+dem2pts(on_offshore25_dem_name,
+        easting_min=project_slide.eastingmin25,
+        easting_max=project_slide.eastingmax25,
+        northing_min=project_slide.northingmin25,
+        northing_max= project_slide.northingmax25,
+        use_cache=True, verbose=True)
+dem2pts(on_offshore100_dem_name,
+        easting_min=project_slide.eastingmin100,
+        easting_max=project_slide.eastingmax100,
+        northing_min=project_slide.northingmin100,
+        northing_max= project_slide.northingmax100,
+        use_cache=True, verbose=True)
+dem2pts(nsw_dem_name,
+        easting_min=project_slide.eastingmin_nsw,
+        easting_max=project_slide.eastingmax_nsw,
+        northing_min=project_slide.northingmin_nsw,
+        northing_max= project_slide.northingmax_nsw,
+        use_cache=True, verbose=True)
+
+print 'create offshore'
+G11 = Geospatial_data(file_name = project_slide.offshore_dem_name1 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name2 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name3 + '.xya')
+G12 = Geospatial_data(file_name = project_slide.offshore_dem_name4 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name5 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name6 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name7 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name8 + '.xya')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name9 + '.xya')
+print 'create onshore'
+G2 = Geospatial_data(file_name = project_slide.on_offshore25_dem_name + '.pts')
+G3 = Geospatial_data(file_name = project_slide.on_offshore100_dem_name + '.pts')
+G4 = Geospatial_data(file_name = project_slide.nsw_dem_name + '.pts')
+print 'add'
+G = G11.clip(Geospatial_data(project_slide.poly_surveyclip)) +\
+    G12.clip(Geospatial_data(project_slide.polyAll)) +\
+    G2.clip(Geospatial_data(project_slide.poly_25mclip)) +\
+    G3.clip(Geospatial_data(project_slide.poly_origsyd)) +\
+    (G4.clip(Geospatial_data(project_slide.polyAll))).clip_outside(Geospatial_data(project_slide.poly_surveyclip)).clip_outside(Geospatial_data(project_slide.poly_origsyd))
+print 'export points'
+G.export_points_file(project_slide.combined_dem_name + '.pts')
+#G.export_points_file(project_slide.combined_dem_name + '.xya')
+
+#----------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with a tagged
+# boundary and interior regions defined in project_slide.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+remainder_res = 1000000.
+local_res = 50000.
+coast_res = 50000.
+interior_regions = [[project_slide.poly_syd1, local_res],
+                    [project_slide.poly_coast, coast_res]]
+
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project_slide.polyAll,
+           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
+                              'e3': [3], 'e4':[4]},
+           'maximum_triangle_area': remainder_res,
+           'filename': meshname,
+           'interior_regions': interior_regions},
+          verbose = True, evaluate=False)
+print 'created mesh'
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+domain = Domain(meshname, use_cache = True, verbose = True)
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+print domain.statistics()
+
+domain.set_name(project_slide.basename)
+domain.set_datadir(project_slide.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_slide.combined_dem_name + '.pts',
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.1
+                    )
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+
+    
+from smf import slide_tsunami    
+tsunami_source = slide_tsunami(length=project_slide.shovel_length,
+                               width=project_slide.shovel_width,
+                               depth=project_slide.shovel_depth,
+                               slope=project_slide.shovel_slope,
+                               thickness=project_slide.shovel_thickness, 
+                               x0=project_slide.slide_origin_shovel_b[0], 
+                               y0=project_slide.slide_origin_shovel_b[1], 
+                               alpha=project_slide.shovel_alpha, 
+                               domain=domain)
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+
+domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )
+
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+import time
+from Numeric import allclose
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'e2')
+    stagestep = domain.get_quantity('stage') 
+
+    if allclose(t, 30):
+        slide = Quantity(domain)
+        slide.set_values(tsunami_source)
+        domain.set_quantity('stage', slide + stagestep)
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
+print 'finished'
