Index: /anuga_work/production/perth_2006/build_perth.py
===================================================================
--- /anuga_work/production/perth_2006/build_perth.py	(revision 4077)
+++ /anuga_work/production/perth_2006/build_perth.py	(revision 4077)
@@ -0,0 +1,161 @@
+"""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_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.abstract_2d_finite_volumes.util 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
+#-------------------------------------------------------------------------------
+
+
+# topography directory filenames
+onshore_dir_name = project.onshore_dir_name
+coast_dir_name = project.coast_dir_name
+islands_dir_name = project.islands_dir_name
+islands_dir_name1 = project.islands_dir_name1
+islands_dir_name1 = project.islands_dir_name2
+offshore_dir_name = project.offshore_dir_name
+
+# creates DEM from asc data
+convert_dem_from_ascii2netcdf(onshore_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(islands_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(islands_dir_name1, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(islands_dir_name2, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+dem2pts(onshore_dir_name,
+#        easting_min=project.eastingmin,
+#        easting_max=project.eastingmax,
+#        northing_min=project.northingmin,
+#        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
+
+#creates pts file for islands DEM
+dem2pts(islands_dir_name, use_cache=True, verbose=True)
+dem2pts(islands_dir_name1, use_cache=True, verbose=True)
+dem2pts(islands_dir_name2, use_cache=True, verbose=True)
+
+print'create Geospatial data objects from topographies'
+G1 = Geospatial_data(file_name = project.onshore_dir_name + '.pts')
+G2 = Geospatial_data(file_name = project.coast_dir_name + '.xya')
+G3 = Geospatial_data(file_name = project.islands_dir_name + '.pts')
+G4 = Geospatial_data(file_name = project.islands_dir_name1 + '.pts')
+G5 = Geospatial_data(file_name = project.islands_dir_name2 + '.pts')
+G_off = Geospatial_data(file_name = project.offshore_dir_name + '.xya')
+
+
+
+
+
+print 'G_offshore_data_old_e'#, G_offshore_old_e.get_data_points() 
+print'add all geospatial objects'
+G = G1 + G2 + G3 + G4 + G5 + G_off 
+
+#print'clip combined geospatial object by bounding polygon'
+G_clipped = G.clip(project.bounding_polygon)
+#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)
+#G_clipped.export_points_file(project.combined_dir_name + '.pts')
+G_clipped.export_points_file(project.combined_dir_name + '.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/perth_2006/project.py
===================================================================
--- /anuga_work/production/perth_2006/project.py	(revision 4076)
+++ /anuga_work/production/perth_2006/project.py	(revision 4077)
@@ -10,12 +10,10 @@
 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
 
-if sys.platform == 'win32':
-    home = getenv('INUNDATIONHOME')
-    user = getenv('USERPROFILE')
+# file and system info
+#---------------------------------
+codename = 'project.py'
 
-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
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+user = get_user_name()
 
 # INUNDATIONHOME is the inundation directory, not the data directory.
@@ -24,117 +22,112 @@
 #Making assumptions about the location of scenario data
 state = 'western_australia'
+scenario_name = 'perth'
 scenario_dir_name = 'perth_tsunami_scenario_2006'
 
 # onshore data provided by WA DLI
-onshore_name = '' # original
+onshore_name = 'perth_dli_ext' # original
+#island
+island_name = 'rott_dli_ext' # original
+island_name1 = 'gard_dli_ext'
+island_name2 = 'carnac_island_dted'
 
 # AHO + DPI data
-offshore_name1 = 'XY100011610'
-offshore_name2 = 'XY100011611'
-offshore_name3 = 'XY100011613'
-offshore_name4 = 'XY100011614'
-offshore_name5 = 'XY100011616'
-offshore_name6 = 'XY100011617'
-offshore_name7 = 'XY100011618'
-offshore_name8 = 'XY100011621'
-offshore_name9 = 'XY100011623'
-offshore_name10 = 'XY100011745'
-offshore_name11 = 'XY100011746'
-offshore_name12 = 'XY100017530'
-offshore_name13 = 'XY100017532'
-offshore_name14 = 'XY100017538'
-offshore_name15 = 'XY100017540'
-offshore_name16 = 'XYBR66'
-offshore_name17 = 'XYBR70'
-offshore_name18 = 'XYBR80'
-offshore_name19 = 'XYBR88'
-offshore_name20 = 'XYBR93'
-offshore_name21 = 'XYBR0110'
-offshore_name22 = 'XYWADPI'
+coast_name = 'waterline'
+offshore_name = 'perth_bathymetry'
 
-# developed by NM&I
-coast_name = 'coastline'
+#final topo name
+combined_name ='perth_combined_elevation'
 
-boundary_basename = 'SU-AU' # Mw ?
+topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep
+topographies_time_dir = topographies_dir+build_time+sep
 
-#swollen/ all data output
-basename = 'source'
-codename = 'project.py'
+#input topo file location
+onshore_dir_name = topographies_dir + onshore_name
+island_dir_name = topographies_dir + island_name
+island_dir_name1 = topographies_dir + island_name1
+island_dir_name2 = topographies_dir + island_name2
+
+coast_dir_name = topographies_dir + coast_name
+offshore_dir_name = topographies_dir + offshore_name
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_time_dir_name = topographies_time_dir + combined_name
+combined_time_dir_final_name = topographies_time_dir + combined_final_name
+
 
 #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 + 'perth_gauges.csv'
-community_filename = gaugedir + 'CHINS_v2.csv'
-community_perth = gaugedir + 'community_perth.csv'
-codedir = getcwd()+sep                            
-codedirname = codedir + 'project.py'
-meshname = outputtimedir + 'mesh_' + basename
+meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
 
-# Necessary if using point datasets, rather than grid
-onshore_dem_name = datadir + onshore_name
-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
-offshore_dem_name10 = datadir + offshore_name10
-offshore_dem_name11 = datadir + offshore_name11
-offshore_dem_name12 = datadir + offshore_name12
-offshore_dem_name13 = datadir + offshore_name13
-offshore_dem_name14 = datadir + offshore_name14
-offshore_dem_name15 = datadir + offshore_name15
-offshore_dem_name16 = datadir + offshore_name16
-offshore_dem_name17 = datadir + offshore_name17
-offshore_dem_name18 = datadir + offshore_name18
-offshore_dem_name19 = datadir + offshore_name19
-offshore_dem_name20 = datadir + offshore_name20
-offshore_dem_name21 = datadir + offshore_name21
-offshore_dem_name22 = datadir + offshore_name22
+polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep
+tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep
 
-coast_dem_name = datadir + coast_name
+#boundaries locations
+boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
+boundaries_in_dir_name = boundaries_in_dir + boundaries_name
+boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + boundaries_name
+#boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
+#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
+#ideas
+#boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep
 
-combined_dem_name   = datadir + 'perth_combined_elevation'
+#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
+
+#output locations
+output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
+output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep
+output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+gauge_name = 'perth.csv'
+gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+
+
+
+
+
+
+
 
 ###############################
 # Domain definitions
 ###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
 
-# bounding box for clipping MOST/URS output (much bigger than study area)
-##south = degminsec2decimal_degrees(-19,0,0)
-##north = degminsec2decimal_degrees(-17,15,0)
-##west  = degminsec2decimal_degrees(120,0,0)
-##east  = degminsec2decimal_degrees(122,30,0) 
-##
-##d0 = [south, west]
-##d1 = [south, east]
-##d2 = [north, east]
-##d3 = [north, west]
-##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
-##refzone = zone
+bounding_polygon = read_polygon(polygons_dir+'bounding_poly.csv')
 
-# bounding polygon for study area
-polyAll = read_polygon(polygondir+'extent.csv')
+refzone = 50 
 
-# plot bounding polygon and make sure BC info surrounds it
-#plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
-print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
+print "bounding_polygon", bounding_polygon
+print 'poly area', polygon_area(bounding_polygon)/1000000.0
+
+
+#Interior regions
+
+#cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+#assert zone == refzone
+
+poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20pts.csv')
+
+poly_cbd = read_polygon(polygons_dir+'cbd_pts.csv')
+
+poly_penguin = read_polygon(polygons_dir+'penguin_pts.csv')
+assert zone == refzone
+
 
 ###################################################################
 # Clipping regions for export to asc and regions for clipping data
 ###################################################################
-
+'''
 # exporting asc grid
 eastingmin = 406215.87
@@ -158,20 +151,4 @@
 print 'Area of coastal polygon', polygon_area(poly_perth3)/1000000.0
 #print 'Area of cable beach polygon', polygon_area(poly_perth4)/1000000.0
+'''
 
-for i in poly_perth3:
-    v = is_inside_polygon(i,poly_perth1, verbose=False)
-    if v == False: print v
-
-def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
-    from anuga.utilities.polygon import polygon_area
-    
-    # TO DO check if any of the regions fall inside one another
-    no_triangles = 0.0
-    area = polygon_area(bounding_poly)
-    for i,j in interior_regions:
-        this_area = polygon_area(i)
-        no_triangles += this_area/j
-        area -= this_area
-        print j, this_area/1000000., area/1000000.
-    no_triangles += area/remainder_res
-    return int(no_triangles/0.7)
