Index: anuga_work/production/dampier_2006/project_urs.py
===================================================================
--- anuga_work/production/dampier_2006/project_urs.py	(revision 4542)
+++ anuga_work/production/dampier_2006/project_urs.py	(revision 4542)
@@ -0,0 +1,294 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd,umask
+from os.path import expanduser, basename
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
+import sys
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
+from time import localtime, strftime, gmtime
+from anuga.utilities.system_tools import get_user_name, get_host_name
+
+# file and system info
+#---------------------------------
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir   
+user = get_user_name()
+host = get_host_name()
+#needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
+umask(002)
+
+#time stuff
+time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
+build_time = time+'_build'
+run_time = time+'_run'
+
+#tide = -2.5
+#tide = 0.0
+tide = 2.4
+
+#Maybe will try to make project a class to allow these parameters to be passed in.
+alpha = 0.1
+friction=0.01
+starttime=3600
+finaltime=25000
+setup='final'
+source='onslow'
+
+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=1
+    yieldstep=60
+if setup =='store':
+    print'store'
+    res_factor=10
+    time_thinning=48
+    yieldstep=240    
+
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
+
+
+#Making assumptions about the location of scenario data
+state = 'western_australia'
+scenario_name = 'dampier'
+scenario = 'dampier_tsunami_scenario_2006'
+
+# onshore data provided by WA DLI
+onshore_name = 'DLI_DTED_raster_clipped' # original
+
+# AHO + DPI data + colin French coastline
+coast_name = 'coastline_edited_w_DEM'
+offshore_name = 'clipped_bathy'
+offshore1_name = 'elev_501'
+offshore2_name = 'inferrec_e'
+
+#final topo name
+combined_name ='dampier_combined_elevation'
+combined_smaller_name = 'dampier_combined_elevation_small'
+combined_smallest_name = 'dampier_combined_elevation_smallest'
+
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
+
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'070112'+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
+coast_in_dir_name = topographies_in_dir + coast_name
+offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore1_in_dir_name = topographies_in_dir + offshore1_name
+offshore2_in_dir_name = topographies_in_dir + offshore2_name
+
+onshore_dir_name = topographies_dir + onshore_name
+coast_dir_name = topographies_dir + coast_name
+offshore_dir_name = topographies_dir + offshore_name
+offshore1_dir_name = topographies_dir + offshore1_name
+offshore2_dir_name = topographies_dir + offshore2_name
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_dir_name = topographies_dir + combined_smaller_name
+combined_smallest_dir_name = topographies_dir + combined_smallest_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+'2007polys'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+
+if source =='dampier':
+    boundaries_name = 'pt_hedland_3854_17042007' #Dampier gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep
+
+if source=='onslow':
+    boundaries_name = 'pt_hedland_3859_16052007' #onslow_hedland_broome gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
+    
+if source=='exmouth':
+    boundaries_name = 'pt_hedland_3103_18052007' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
+
+#boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'1_10000'+sep
+boundaries_in_dir_name = boundaries_in_dir + boundaries_name
+boundaries_dir = anuga_dir+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + boundaries_name
+
+#output locations
+output_dir = anuga_dir+'outputs'+sep
+output_build_time_dir = output_dir+build_time+'_'+source+sep
+output_run_time_dir = output_dir +run_time+dir_comment+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+gauge_name = 'dampier_gauges_up2.csv' #'dampier.csv'
+gauge_name_simple = 'dampier_gauges_simple.csv' 
+gauge_name_test = 'dampier_gauges_up_test.csv' #'dampier.csv'
+gauges_dir = anuga_dir+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+gauges_dir_name_simple = gauges_dir + gauge_name_simple
+gauges_dir_name_test = gauges_dir + gauge_name_test
+
+community_filename = gauges_dir + 'CHINS_v2.csv'
+community_KD = gauges_dir + 'community_KD.csv'
+
+buildings_filename = gauges_dir + 'Dampier_res_Project.csv'
+buildings_filename_out = 'Dampier_res_Project_modified.csv'
+###############################
+# Domain definitions
+###############################
+
+#refzone = 50 
+south_boundary = degminsec2decimal_degrees(-20,58,0)
+#north_boundary = degminsec2decimal_degrees(-20,13,0)
+north_boundary = degminsec2decimal_degrees(-19,30,0)
+west_boundary = degminsec2decimal_degrees(116,15,0)
+east_boundary = degminsec2decimal_degrees(117,11,0)
+
+poly_all = read_polygon(polygons_dir+'extent.csv')
+print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0
+
+res_poly_all = 150000*res_factor
+#res_poly_all = 150000
+
+poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')
+res_pipeline = 2000
+
+###############################
+# Interior region definitions
+###############################
+
+poly_region = read_polygon(polygons_dir+'region.csv')
+res_region = 50000*res_factor
+
+poly_dampier = read_polygon(polygons_dir+'dampier_town.csv')
+res_dampier = 500*res_factor
+
+poly_karratha = read_polygon(polygons_dir+'karrathav2.csv')
+res_karratha = 15000*res_factor
+
+poly_karratha_town = read_polygon(polygons_dir+'karratha_townv2.csv')
+res_karratha_town = 500*res_factor
+
+poly_facility = read_polygon(polygons_dir+'facility.csv')
+res_facility = 1000*res_factor
+
+poly_delambre = read_polygon(polygons_dir+'delambre.csv')
+res_delambre = 1000*res_factor
+
+poly_coast = read_polygon(polygons_dir+'coastpoly.csv')
+res_coast = 5000*res_factor
+
+poly_NWislands = read_polygon(polygons_dir+'nw_islands_area.csv')
+res_NWislands = 50000*res_factor
+
+poly_island0 = read_polygon(polygons_dir+'island0.csv')
+res_island0 = res_poly_all
+
+poly_island1 = read_polygon(polygons_dir+'island1.csv')
+res_island0 = res_poly_all
+
+poly_island2 = read_polygon(polygons_dir+'island2.csv')
+res_island0 = res_poly_all
+
+poly_island3 = read_polygon(polygons_dir+'island3.csv')
+res_island0 = res_poly_all
+
+res_islands = 5000*res_factor
+
+poly_ref_nw4 = read_polygon(polygons_dir+'ref_nw4.csv')
+res_ref_nw4 = res_islands
+
+poly_island4 = read_polygon(polygons_dir+'island4.csv')
+res_island0 = res_poly_all
+
+poly_ref_nw5 = read_polygon(polygons_dir+'ref_nw5.csv')
+res_ref_nw5 = res_islands
+
+poly_island5 = read_polygon(polygons_dir+'island5.csv')
+res_island0 = res_poly_all
+
+poly_ref_nw6 = read_polygon(polygons_dir+'ref_nw6.csv')
+res_ref_nw6 = res_islands
+
+poly_island6 = read_polygon(polygons_dir+'island6.csv')
+res_island0 = res_poly_all
+
+poly_ref_nw7 = read_polygon(polygons_dir+'ref_nw7.csv')
+res_ref_nw7 = res_islands
+
+poly_island7 = read_polygon(polygons_dir+'island7.csv')
+res_island0 = res_poly_all
+
+poly_ref_nw8 = read_polygon(polygons_dir+'ref_nw8.csv')
+res_ref_nw8 = res_islands
+
+poly_island8 = read_polygon(polygons_dir+'island8.csv')
+res_island0 = res_poly_all
+
+
+#plot_polygons([poly_dampier,poly_karratha,poly_karratha_town,poly_delambre,
+#                poly_coast,poly_NWislands,poly_island0,poly_island1,poly_island2,
+#                poly_island3,poly_island4,poly_island5,poly_island6,
+#                poly_island7,poly_island8,poly_ref_nw4,poly_ref_nw5,
+#                poly_ref_nw6,poly_ref_nw7,poly_ref_nw8,poly_all],'poly_pic')
+
+interior_regions = [[poly_dampier,res_dampier], 
+                    [poly_karratha,res_karratha],[poly_karratha_town,res_karratha_town],
+                    [poly_delambre,res_delambre],[poly_coast,res_coast],
+                    [poly_facility,res_facility],
+                    #[poly_NWislands,res_NWislands],
+                    [poly_island0,res_island0],[poly_island1,res_island0],
+                    [poly_island2,res_island0],[poly_island3,res_island0],
+                    [poly_island4,res_island0],[poly_island5,res_island0],
+                    [poly_island6,res_island0],[poly_island7,res_island0],
+                    [poly_island8,res_island0],[poly_ref_nw4,res_ref_nw4],
+                    [poly_ref_nw5,res_ref_nw5],[poly_ref_nw6,res_ref_nw6],
+                    [poly_ref_nw7,res_ref_nw7],[poly_ref_nw8,res_ref_nw8]
+#                    ,[poly_pipeline,res_pipeline]
+                    ]
+
+interior_regions_test = [[poly_dampier,res_dampier], 
+                    [poly_karratha,res_karratha],[poly_karratha_town,res_karratha_town]]
+
+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
+###################################################################
+
+poly_mainland = read_polygon(polygons_dir+'mainland_only.csv')
+
+# exporting asc grid - Dampier gas facility
+e_min_area = 474000
+e_max_area = 480000
+n_min_area = 7719000
+n_max_area = 7725000
+
+# residential Dampier
+e_min_area_d = 468450
+e_max_area_d = 470820
+n_min_area_d = 7714050
+n_max_area_d = 7716625
+
+# residential Karratha
+e_min_area_k = 474690
+e_max_area_k = 489900
+n_min_area_k = 7704760
+n_max_area_k = 7712600
Index: anuga_work/production/dampier_2006/run_dampier_urs.py
===================================================================
--- anuga_work/production/dampier_2006/run_dampier_urs.py	(revision 4542)
+++ anuga_work/production/dampier_2006/run_dampier_urs.py	(revision 4542)
@@ -0,0 +1,261 @@
+"""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_urs.py
+The output sww file is stored in project_urs.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project_urs.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,umask
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Field_boundary
+from Numeric import allclose
+from anuga.shallow_water.data_manager import export_grid
+
+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
+# Application specific imports
+import project_urs                 # Definition of file names and polygons
+
+def run_model(**kwargs):
+    
+    tide = kwargs['tide']
+    alpha = kwargs['alpha']
+    friction = kwargs['friction']
+    time_thinning = kwargs['time_thinning']
+    scenario_name = kwargs['aa_scenario_name']
+    
+    #------------------------------------------------------------------------------
+    # Copy scripts to time stamped output directory and capture screen
+    # output to file
+    #------------------------------------------------------------------------------
+
+    #copy script must be before screen_catcher
+    print 'tide',tide
+    kwargs['est_num_trigs']=project_urs.trigs_min
+    kwargs['num_cpu']=numprocs
+    kwargs['host']=project_urs.host
+    kwargs['res_factor']=project_urs.res_factor
+    kwargs['starttime']=project_urs.starttime
+    kwargs['yieldstep']=project_urs.yieldstep
+    kwargs['finaltime']=project_urs.finaltime
+    
+    kwargs['output_dir']=project_urs.output_run_time_dir
+    kwargs['bathy_file']=project_urs.combined_dir_name + '.txt'
+#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
+#    kwargs['Completed']=''
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print 'output_dir',kwargs['output_dir']
+    if myid == 0:
+        copy_code_files(kwargs['output_dir'],__file__, 
+                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
+
+        store_parameters(**kwargs)
+
+    barrier()
+
+    print "Processor Name:",get_processor_name()
+
+    # creates copy of code in output dir
+    print 'min triangles', project_urs.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_urs.py along with
+    # resolutions (maximal area of per triangle) for each polygon
+    #--------------------------------------------------------------------------
+
+    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
+    # causes problems with the ability to cache set quantity which takes alot of times
+    if myid == 0:
+    
+        print 'start create mesh from regions'
+
+        create_mesh_from_regions(project_urs.poly_all,
+                             boundary_tags={'back': [2,3], 'side': [0, 1, 4],
+                                            'ocean': [5]},
+                             maximum_triangle_area=project_urs.res_poly_all,
+                             interior_regions=project_urs.interior_regions,
+                             filename=project_urs.meshes_dir_name+'.msh',
+                             use_cache=False,
+                             verbose=True)
+    barrier()
+
+    #-------------------------------------------------------------------------
+    # 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_urs.meshes_dir_name+'.msh', use_cache=False, verbose=True)
+       
+    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_urs.poly_mainland, -1.0)], default = tide,
+                                 geo_reference = domain.geo_reference)
+        domain.set_quantity('stage', IC)
+        domain.set_quantity('friction', friction) 
+        
+        print 'Start Set quantity'
+
+        domain.set_quantity('elevation', 
+                            filename = kwargs['bathy_file'],
+                            use_cache = True,
+                            verbose = True,
+                            alpha = 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(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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+    #domain.limit2007 = 1
+
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
+
+    Bf = Field_boundary(kwargs['boundary_file'],
+                    domain, time_thinning=time_thinning, mean_stage=tide, 
+                    use_cache=True, verbose=True)
+                    
+    kwargs['input_start_time']=domain.starttime
+
+    print 'finished reading boundary file'
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([tide,0,0])
+
+    print'set_boundary'
+
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': Bf}) 
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #----------------------------------------------------------------------------
+
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+
+    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
+                       ,skip_initial_step = True): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+    
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = True): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')   
+
+    x, y = domain.get_maximum_inundation_location()
+    q = domain.get_maximum_inundation_elevation()
+
+    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+
+    print 'That took %.2f seconds' %(time.time()-t0)
+    
+    kwargs['completed']=str(time.time()-t0)
+    
+    if myid == 0:
+        store_parameters(**kwargs)
+
+    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
+    export_grid(swwfile, extra_name_out = 'town',
+                quantities = ['elevation','depth','stage','speed'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+                timestep = None,
+                reduction = max,
+                cellsize = 25,
+                NODATA_value = -9999,
+                easting_min = project_urs.eastingmin,
+                easting_max = project_urs.eastingmax,
+                northing_min = project_urs.northingmin,
+                northing_max = project_urs.northingmax,
+                verbose = True,
+                origin = None,
+                datum = 'WGS84',
+                format = 'asc')
+    
+    buildings_filename = project_urs.buildings_filename
+    buildings_filename_out = project_urs.buildings_filename_out
+                
+    inundation_damage(swwfile+'.sww', buildings_filename, buildings_filename_out)
+    print '\n Augmented building file written to %s \n' %buildings_filename_out
+    
+    barrier()
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+    run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name,
+              ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user,
+              alpha = project_urs.alpha, friction=project_urs.friction,
+              time_thinning = project_urs.time_thinning,
+              dir_comment=project_urs.dir_comment)
+
+
