Index: /anuga_work/production/carnarvon/build_boundary.py
===================================================================
--- /anuga_work/production/carnarvon/build_boundary.py	(revision 5755)
+++ /anuga_work/production/carnarvon/build_boundary.py	(revision 5755)
@@ -0,0 +1,99 @@
+from anuga.utilities.numerical_tools import ensure_numeric
+from Scientific.IO.NetCDF import NetCDFFile
+from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
+    compress,zeros,fabs,allclose,ones
+from anuga.utilities.polygon import inside_polygon,read_polygon
+from os import sep
+from time import localtime, strftime, gmtime
+from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
+import os
+import project
+
+######################################
+# Create urs boundary from mux2files #
+######################################
+
+dir=os.path.join(project.muxhome,'mux')
+
+# taken from David's event list
+urs_filenames = [
+    os.path.join(dir,'Java-0000-z.grd'),
+    os.path.join(dir,'Java-0001-z.grd'),
+    os.path.join(dir,'Java-0002-z.grd'),
+    os.path.join(dir,'Java-0003-z.grd'),
+    os.path.join(dir,'Java-0005-z.grd'), 
+    os.path.join(dir,'Java-0006-z.grd'), 
+    os.path.join(dir,'Java-0007-z.grd'), 
+    os.path.join(dir,'Java-0008-z.grd'), 
+    os.path.join(dir,'Java-0010-z.grd'), 
+    os.path.join(dir,'Java-0011-z.grd'), 
+    os.path.join(dir,'Java-0012-z.grd'), 
+    os.path.join(dir,'Java-0013-z.grd'),  
+    os.path.join(dir,'Java-0015-z.grd'),  
+    os.path.join(dir,'Java-0016-z.grd'),  
+    os.path.join(dir,'Java-0017-z.grd'),  
+    os.path.join(dir,'Java-0018-z.grd'),  
+    os.path.join(dir,'Java-0020-z.grd'),  
+    os.path.join(dir,'Java-0021-z.grd'),  
+    os.path.join(dir,'Java-0022-z.grd'),  
+    os.path.join(dir,'Java-0023-z.grd'),  
+    os.path.join(dir,'Java-0025-z.grd'),  
+    os.path.join(dir,'Java-0026-z.grd'),  
+    os.path.join(dir,'Java-0027-z.grd'),  
+    os.path.join(dir,'Java-0028-z.grd'),  
+    os.path.join(dir,'Java-0029-z.grd'),  
+    os.path.join(dir,'Java-0030-z.grd'),  
+    os.path.join(dir,'Java-0031-z.grd'),  
+    os.path.join(dir,'Java-0032-z.grd'),  
+    os.path.join(dir,'Java-0033-z.grd'),  
+    os.path.join(dir,'Java-0034-z.grd'),  
+    os.path.join(dir,'Java-0035-z.grd'),  
+    os.path.join(dir,'Java-0036-z.grd'),  
+    os.path.join(dir,'Java-0037-z.grd'),  
+    os.path.join(dir,'Java-0038-z.grd'),  
+    os.path.join(dir,'Java-0039-z.grd'),  
+    os.path.join(dir,'Java-0040-z.grd'),  
+    os.path.join(dir,'Java-0041-z.grd'),  
+    os.path.join(dir,'Java-0042-z.grd'),  
+    os.path.join(dir,'Java-0043-z.grd'),  
+    os.path.join(dir,'Java-0044-z.grd'),  
+    os.path.join(dir,'Java-0045-z.grd'),  
+    os.path.join(dir,'Java-0046-z.grd'),  
+    os.path.join(dir,'Java-0047-z.grd'),
+    os.path.join(dir,'Java-0048-z.grd')]  
+ 
+print 'number of sources', len(urs_filenames)
+
+#from os import sys
+#sys.exit()
+# for all mux files
+#urs_filenames={}
+##for filename in os.listdir(dir)
+##    basename = filename[:-7]
+##
+##    urs_filenames[basename] = True
+    
+# Make list of basenames  
+#urs_filenames = [os.path.join(dir, basename) for basename in urs_filenames.keys()]
+
+# AS per David Burbidge email on Friday 4th July the mag 9.3 event
+# has 1m worth of slip on each sub fault therefore mutliple each unit
+# source by the slip (10.4544) and sum the 44 time series together to
+# get the time series for this event at the points on your boundary. 
+
+weight_factor = 10.4544
+weights=weight_factor*ones(len(urs_filenames),Float)
+
+scenario_name=project.scenario_name
+order_filename=os.path.join(project.boundaries_dir, 'thinned_boundary_ordering.txt')
+
+print 'reading', order_filename
+# Create ordered sts file
+print 'creating sts file'
+
+urs2sts(urs_filenames,
+        basename_out=os.path.join(project.boundaries_dir,scenario_name),
+        ordering_filename=order_filename,
+        weights=weights,
+        mean_stage=project.tide,
+        verbose=True)
Index: /anuga_work/production/carnarvon/build_carnarvon.py
===================================================================
--- /anuga_work/production/carnarvon/build_carnarvon.py	(revision 5754)
+++ /anuga_work/production/carnarvon/build_carnarvon.py	(revision 5755)
@@ -27,5 +27,5 @@
 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
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
 
 # Application specific imports
@@ -37,14 +37,10 @@
 #------------------------------------------------------------------------------
 
-start_screen_catcher(project.output_build_time_dir)
-
-print 'time stamp: ',project.gtime
-print 'USER: ', project.user
-
-
 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
 
 #-------------------------------------------------------------------------------
@@ -55,20 +51,25 @@
 # Fine pts file to be clipped to area of interest
 #-------------------------------------------------------------------------------
+print"project.poly_all",project.poly_all
+print"project.combined_dir_name",project.combined_dir_name
 
 # topography directory filenames
 onshore_in_dir_name = project.onshore_in_dir_name
 coast_in_dir_name = project.coast_in_dir_name
-#island_in_dir_name = project.island_in_dir_name
 offshore_in_dir_name = project.offshore_in_dir_name
+offshore_in_dir_name1 = project.offshore_in_dir_name1
+offshore_in_dir_name2 = project.offshore_in_dir_name2
 
 onshore_dir_name = project.onshore_dir_name
 coast_dir_name = project.coast_dir_name
-#island_dir_name = project.island_dir_name
+
 offshore_dir_name = project.offshore_dir_name
+offshore_dir_name1 = project.offshore_dir_name1
+offshore_dir_name2 = project.offshore_dir_name2
+
 
 # creates DEM from asc data
 print "creates DEMs from ascii data"
 convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
-#convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
 
 #creates pts file for onshore DEM
@@ -86,78 +87,36 @@
 
 print'create Geospatial data1 objects from topographies'
-G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
-G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
-#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
-G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True)
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True)
 
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
+G_off1 = Geospatial_data(file_name = offshore_in_dir_name1 + '.txt',verbose=True)
+G_off2 = Geospatial_data(file_name = offshore_in_dir_name2 + '.txt',verbose=True)
 print'add all geospatial objects'
-G = G1 + G2 + G_off 
+G = G1 + G2 + G_off + G_off1 + G_off2
 
 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
+G_clipped = G.clip(project.poly_all)
 
 print'export combined DEM file'
 if access(project.topographies_dir,F_OK) == 0:
     mkdir (project.topographies_dir)
-G.export_points_file(project.combined_dir_name + '.xya')
-#G_clipped.export_points_file(project.combined_dir_name + '.xya')
+G_clipped.export_points_file(project.combined_dir_name + '.txt')
 
-'''
-print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
-G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
-print'split'
-G_all_1, G_all_2 = G_all.split(.10)
-print'export 1'
-G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
-print'export 2'
-G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
+print'project.combined_dir_name + .txt',project.combined_dir_name + '.txt'
 
-
-#-------------------------------------------------------------------------
-# 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')
-
-'''
-
-
-
-
-
-
-
+###-------------------------------------------------------------------------
+### Convert URS to SWW file for boundary conditions 
+###-------------------------------------------------------------------------
+##print 'starting to create boundary conditions'
+##from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
+##
+##boundaries_in_dir_name = project.boundaries_in_dir_name
+##print 'boundaries_in_dir_name',project.boundaries_in_dir_name
+##
+##
+###import sys; sys.exit()
+##
+##urs_ungridded2sww(project.boundaries_in_dir_name, project.boundaries_in_dir_name,
+##                  verbose=True, mint=4000, maxt=80000, zscale=1)
+##
Index: /anuga_work/production/carnarvon/export_results.py
===================================================================
--- /anuga_work/production/carnarvon/export_results.py	(revision 5755)
+++ /anuga_work/production/carnarvon/export_results.py	(revision 5755)
@@ -0,0 +1,81 @@
+import project, os
+import sys
+
+from anuga.shallow_water.data_manager import sww2dem
+from os import sep
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+
+
+time_dir = '20080908_154507_run_final_0.0_polyline_alpha0.1_kvanputt'
+
+
+cellsize = 25
+#cellsize = 150
+timestep = 0
+directory = project.output_dir
+name = directory+sep+time_dir+sep+project.scenario_name
+
+
+
+
+is_parallel = False
+#is_parallel = True
+
+if is_parallel == True: nodes = 4
+print 'output dir:', name
+
+
+var = [0,4] 
+#var = [2,3] # depth and Speed
+#var = [2,3] # elevation, depth and Speed
+
+
+for which_var in var:
+    if which_var == 0:  # Stage
+        outname = name + '_stage'
+        quantityname = 'stage'
+
+    if which_var == 1:  # Absolute Momentum
+        outname = name + '_momentum'
+        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'
+        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:
+    #    print 'is_parallel',is_parallel
+        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 = timestep,
+                        cellsize = cellsize,      
+                        easting_min = project_grad.e_min_area,
+                        easting_max = project_grad.e_max_area,
+                        northing_min = project_grad.n_min_area,
+                        northing_max = project_grad.n_max_area,        
+                        reduction = max, 
+                        verbose = True,
+                        format = 'asc')
+    else:
+        print 'start sww2dem'
+        sww2dem(name, basename_out = outname,
+                    quantity = quantityname,
+                    timestep = timestep,
+                    cellsize = cellsize,
+                    number_of_decimal_places = 4,
+                    reduction = max, 
+                    verbose = True,
+                    format = 'asc')
+
Index: /anuga_work/production/carnarvon/get_gauges.py
===================================================================
--- /anuga_work/production/carnarvon/get_gauges.py	(revision 5755)
+++ /anuga_work/production/carnarvon/get_gauges.py	(revision 5755)
@@ -0,0 +1,60 @@
+import os
+from os import sep
+import project
+try:
+    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
+except:
+    print 'Cannot import pylab plotting will not work. Csv files are still created'
+from anuga.utilities.polygon import inside_polygon
+from anuga.utilities.numerical_tools import ensure_numeric
+from Scientific.IO.NetCDF import NetCDFFile
+
+def get_sts_gauge_data(filename,verbose=False):
+    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
+        compress,zeros,fabs,take
+    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
+    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
+    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
+    points=transpose(asarray([x.tolist(),y.tolist()]))
+    time=fid.variables['time'][:]+fid.starttime
+    elevation=fid.variables['elevation'][:]
+    print 'hello', fid.variables.keys()
+
+    
+    basename='sts_gauge'
+    quantity_names=['stage','xmomentum','ymomentum']
+    quantities = {}
+    for i, name in enumerate(quantity_names):
+        quantities[name] = fid.variables[name][:]
+
+    maxname = 'max_sts_stage.csv'
+    fid_max = open(project.boundaries_dir+sep+maxname,'w')
+    s = 'x, y, max_stage \n'
+    fid_max.write(s)    
+    for j in range(len(x)):
+        locx=int(x[j])
+        locy=int(y[j])
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
+        fid_max.write(s)
+        
+        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
+        s = 'time, stage, xmomentum, ymomentum \n'
+        fid_sts.write(s)
+        
+        for k in range(len(time)-1):
+            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
+            fid_sts.write(s)
+
+        fid_sts.close()      
+
+    fid.close()
+        
+    return quantities,elevation,time
+
+quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)
+
+print len(elevation), len(quantities['stage'][0,:])
Index: /anuga_work/production/carnarvon/project.py
===================================================================
--- /anuga_work/production/carnarvon/project.py	(revision 5754)
+++ /anuga_work/production/carnarvon/project.py	(revision 5755)
@@ -3,5 +3,5 @@
 """
 
-from os import sep, environ, getenv, getcwd ,umask
+from os import sep, environ, getenv, getcwd
 from os.path import expanduser
 import sys
@@ -10,22 +10,24 @@
 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
 from anuga.utilities.system_tools import get_user_name, get_host_name
+from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
 
 # file and system info
 #---------------------------------
+#codename = 'project.py'
 
-home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir   
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
+muxhome = getenv('MUXHOME')
 user = get_user_name()
 host = get_host_name()
+
 # INUNDATIONHOME is the inundation directory, not the data directory.
-
-#needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
-umask(002)
 
 #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'
-
-tide = 0.6
+print 'gtime: ', gtime
 
 #Making assumptions about the location of scenario data
@@ -34,14 +36,14 @@
 scenario = 'carnarvon_tsunami_scenario'
 
-#Maybe will try to make project a class to allow these parameters to be passed in.
+
+tide = 0.0 #1.0
+
 alpha = 0.1
 friction=0.01
-starttime=10000
-midtime=21600
-finaltime=25000
-export_cellsize=50
+starttime=0
+finaltime=1000 
+export_cellsize=25
 setup='final'
-source='other'
-
+source='polyline'
 
 if setup =='trial':
@@ -61,39 +63,45 @@
     yieldstep=60
 
-dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
+
+# onshore data provided by WA DLI
+onshore_name = 'bathy250_clipland' # original
+
+# AHO + DPI data + colin French coastline
+coast_name = 'DPI_coastlineP'
+offshore_name = 'Shark_Bay_Clip'
+offshore_name1 = 'XYAHD_Clip'
+offshore_name2 = 'DPI_data'
 
 
-# onshore data provided by WA DLI
-onshore_name = 'DLI_orthophoto_DEM' # original
-#islands
-
-# AHO + DPI data
-coast_name = '100k_coastline'
-offshore_name = 'Bathymetry'
-
 #final topo name
-combined_name ='carnarvon_combined_elevation.pts'
-combined_smaller_name = 'carnarvon_combined_elevation_smaller'
+combined_name = scenario_name+'_combined_elevation'
+combined_smaller_name = scenario_name+'_combined_elevation_smaller'
 
 anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
 
-topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
-topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
+topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = anuga_dir+'topographies'+sep
 
-#input topo file location
+# input topo file location
 onshore_in_dir_name = topographies_in_dir + onshore_name
-#island_in_dir_name = topographies_in_dir + island_name
 
 coast_in_dir_name = topographies_in_dir + coast_name
 offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore_in_dir_name1 = topographies_in_dir + offshore_name1
+offshore_in_dir_name2 = topographies_in_dir + offshore_name2
 
 onshore_dir_name = topographies_dir + onshore_name
-#island_dir_name = topographies_dir + island_name
+
 coast_dir_name = topographies_dir + coast_name
 offshore_dir_name = topographies_dir + offshore_name
+offshore_dir_name1 = topographies_dir + offshore_name1
+offshore_dir_name2 = topographies_dir + offshore_name2
 
 #final topo files
 combined_dir_name = topographies_dir + combined_name
-combined_smaller_dir_name = topographies_dir + combined_smaller_name
+#combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
+#combined_time_dir_final_name = topographies_time_dir + combined_final_name
 
 meshes_dir = anuga_dir+'meshes'+sep
@@ -103,54 +111,74 @@
 tide_dir = anuga_dir+'tide_data'+sep
 
-if source =='dampier':
-    boundaries_name = 'broome_3854_17042007' #Dampier gun
-    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep
 
-if source=='onslow':
-    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
-    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
+#boundaries_source = '1'
     
-if source=='exmouth':
-    boundaries_name = 'broome_3103_18052007' #exmouth gun
-    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
+if source=='polyline':
+    boundaries_name = 'perth_3103_28052008' #polyline gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'polyline'+sep+'1_10000'+sep
 
-if source=='other':
-    boundaries_name = 'other' #exmouth gun
-    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
+if source=='test':
+    boundaries_name = 'other' #polyline
+    boundaries_in_dir = anuga_dir+'boundaries'+sep
 
 
 #boundaries locations
-boundaries_in_dir_name = boundaries_in_dir + boundaries_name
+boundaries_in_dir_name = boundaries_in_dir + boundaries_name 
 boundaries_dir = anuga_dir+'boundaries'+sep
-boundaries_dir_name = boundaries_dir + scenario_name
+boundaries_dir_name = boundaries_dir + scenario_name # what it creates???
+boundaries_dir_mux = muxhome
 
 #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 = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
+output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
 
+vertex_filename = output_run_time_dir + 'mesh_vertex.csv'
+
 #gauges
-gauge_name = '???.csv'
+gauge_name = 'carnarvon.csv'
+
 gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+beach_gauges = gauges_dir + 'beach_gauges.csv'
 gauges_dir_name = gauges_dir + gauge_name
 
-buildings_filename = gauges_dir + 'carnarvon_res_Project.csv'
-buildings_filename_out = 'carnarvon_res_Project_modified.csv'
+##buildings_filename = gauges_dir + 'Perth_resA.csv'
+##buildings_filename_out = 'Perth_res_Project_modified.csv'
 
-community_filename = gauges_dir +''
-community_broome = gauges_dir + ''
+###############################
+# Interior region definitions
+###############################
+
+#Initial bounding polygon for data clipping 
+poly_all = read_polygon(polygons_dir+'poly_all.csv')
+res_poly_all = 100000*res_factor
+
+#Polygon designed
+poly_internal_20 = read_polygon(polygons_dir+'Carnarvon20m.csv')
+res_internal_20 = 25000*res_factor
+
+#Polygon designed to 
+poly_internal_5 = read_polygon(polygons_dir+'Carnarvon5m.csv')
+res_internal_5 = 500*res_factor
+
+#Polygon designed to
+poly_internal_10 = read_polygon(polygons_dir+'Carnarvon10m.csv')
+res_internal_10 = 1500*res_factor
+
+#Polygon designed to incorporate
+poly_island_20 = read_polygon(polygons_dir+'Island20m.csv')
+res_island_20 = 50000*res_factor
 
 
-###############################
-# Domain definitions
-###############################
+interior_regions = [[poly_internal_20,res_internal_20],[poly_internal_5,res_internal_5]
+                     ,[poly_island_20,res_island_20],[poly_internal_10,res_internal_10]]
 
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+    
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+print 'min number triangles', trigs_min
+   
 
-# bounding polygon for study area
-poly_all = read_polygon(polygons_dir+'poly_all.csv')
-res_poly_all = 100000*res_factor
+poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
 
 ###################################################################
@@ -158,32 +186,3 @@
 ###################################################################
 
-# exporting asc grid
-eastingmin = 340000
-eastingmax = 350000
-northingmin = 6273400
-northingmax = 6277700
 
-###############################
-# Interior region definitions
-###############################
-
-#digitized polygons
-poly_carnarvon1 = read_polygon(polygons_dir+'neg20_pos10_polygon.csv')
-res_carnarvon1 = 10000*res_factor
-poly_carnarvon2 = read_polygon(polygons_dir+'neg5_pos5_poly_.csv')
-res_carnarvon2 = 500*res_factor
-
-#plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False)
-
-interior_regions = [[poly_carnarvon1,res_carnarvon1],[poly_carnarvon2,res_carnarvon2]]
-
-boundary_tags={'back': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 14],
-               'side': [10,13], 'ocean': [11, 12]}
-
-trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
-
-poly_mainland=read_polygon(polygons_dir+'initial_condition_south.csv')
-
-print 'min number triangles', trigs_min
-
-
Index: /anuga_work/production/carnarvon/run_carnarvon.py
===================================================================
--- /anuga_work/production/carnarvon/run_carnarvon.py	(revision 5754)
+++ /anuga_work/production/carnarvon/run_carnarvon.py	(revision 5755)
@@ -17,4 +17,5 @@
 # Standard modules
 from os import sep
+import os
 from os.path import dirname, basename
 from os import mkdir, access, F_OK
@@ -30,16 +31,20 @@
 from anuga.shallow_water import Field_boundary
 from Numeric import allclose
-from anuga.shallow_water.data_manager import export_grid
+from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
 
 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_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
 from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
 
 # Application specific imports
 import project                 # Definition of file names and polygons
+numprocs = 1
+myid = 0
 
 def run_model(**kwargs):
@@ -62,5 +67,5 @@
         store_parameters(**kwargs)
 
-    barrier()
+   # barrier()
 
     start_screen_catcher(kwargs['output_dir'], myid, numprocs)
@@ -68,11 +73,24 @@
     print "Processor Name:",get_processor_name()
 
-    # filenames
-#    meshes_dir_name = project.meshes_dir_name+'.msh'
-
-    # creates copy of code in output dir
-    print 'min triangles', project.trigs_min,
-    print 'Note: This is generally about 20% less than the final amount'
-
+    
+    #-----------------------------------------------------------------------
+    # Domain definitions
+    #-----------------------------------------------------------------------
+
+    # Read in boundary from ordered sts file
+    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
+
+    # Reading the landward defined points, this incorporates the original clipping
+    # polygon minus the 100m contour
+    landward_bounding_polygon = read_polygon(project.boundaries_dir+'landward_bounding_polygon.txt')
+
+    # Combine sts polyline with landward points
+    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
+    
+    # counting segments
+    N = len(urs_bounding_polygon)-1
+    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4],'ocean': range(N)}
+
+    
     #--------------------------------------------------------------------------
     # Create the triangular mesh based on overall clipping polygon with a
@@ -88,13 +106,17 @@
         print 'start create mesh from regions'
 
-        create_mesh_from_regions(project.poly_all,
-                             boundary_tags=project.boundary_tags,
+        create_mesh_from_regions(bounding_polygon,
+                             boundary_tags=boundary_tags,
                              maximum_triangle_area=project.res_poly_all,
                              interior_regions=project.interior_regions,
                              filename=project.meshes_dir_name+'.msh',
-                             use_cache=False,
+                             use_cache=True,
                              verbose=True)
-    barrier()
-    
+   # barrier()
+
+##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
+##                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
+##                                mesh_file = project.meshes_dir_name+'.msh')
+##        print 'optimal alpha', covariance_value,alpha       
 
     #-------------------------------------------------------------------------
@@ -103,6 +125,4 @@
     print 'Setup computational domain'
 
-    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
-    #above don't work
     domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
     print 'memory usage before del domain',mem_usage()
@@ -123,20 +143,37 @@
         from polygon import Polygon_function
         #following sets the stage/water to be offcoast only
-        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
+        IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
                                  geo_reference = domain.geo_reference)
         domain.set_quantity('stage', IC)
-#        domain.set_quantity('stage', kwargs['tide'])
+        #domain.set_quantity('stage',kwargs['tide'] )
         domain.set_quantity('friction', kwargs['friction']) 
         
-        print 'Start Set quantity',kwargs['bathy_file']
+        print 'Start Set quantity',kwargs['elevation_file']
 
         domain.set_quantity('elevation', 
-                            filename = kwargs['bathy_file'],
-                            use_cache = True,
+                            filename = kwargs['elevation_file'],
+                            use_cache = False,
                             verbose = True,
                             alpha = kwargs['alpha'])
         print 'Finished Set quantity'
-    barrier()
-
+    #barrier()
+
+##    #------------------------------------------------------
+##    # Create x,y,z file of mesh vertex!!!
+##    #------------------------------------------------------
+##        coord = domain.get_vertex_coordinates()
+##        depth = domain.get_quantity('elevation')
+##        
+##        # Write vertex coordinates to file
+##        filename=project.vertex_filename
+##        fid=open(filename,'w')
+##        fid.write('x (m), y (m), z(m)\n')
+##        for i in range(len(coord)):
+##            pt=coord[i]
+##            x=pt[0]
+##            y=pt[1]
+##            z=depth[i]
+##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
+##
 
     #------------------------------------------------------
@@ -157,7 +194,6 @@
     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)
+    domain.tight_slope_limiters = 1
     print 'domain id', id(domain)
-
 
     #-------------------------------------------------------------------------
@@ -166,25 +202,25 @@
     print 'Available boundary tags', domain.get_boundary_tags()
     print 'domain id', id(domain)
-    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
-
-    if project.source != 'other':
-        Bf = Field_boundary(kwargs['boundary_file'],
-                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
-                    use_cache=True, verbose=True)
-                    
-    kwargs['input_start_time']=domain.starttime
-
-    print 'finished reading boundary file'
+   
+    boundary_urs_out=project.boundaries_dir_name
 
     Br = Reflective_boundary(domain)
     Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
-
-
-    print'set_boundary'
+    
+    print 'Available boundary tags', domain.get_boundary_tags()
+    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
+                   domain, mean_stage= project.tide,
+                   time_thinning=1,
+                   default_boundary=Bd,
+                   use_cache=True,
+                   verbose = True,
+                   boundary_polygon=bounding_polygon)
 
     domain.set_boundary({'back': Br,
                          'side': Bd,
                          'ocean': Bf}) 
+
+    kwargs['input_start_time']=domain.starttime
+
     print'finish set boundary'
 
@@ -194,23 +230,8 @@
     t0 = time.time()
 
-    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
+    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = False): 
         domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')     
-
-        if allclose(t, 120):
-            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
-
-        if allclose(t, 720):
-            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
-
-#    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
-#                       ,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')   
+        domain.write_boundary_statistics(tags = 'ocean')
 
     x, y = domain.get_maximum_inundation_location()
@@ -226,36 +247,8 @@
     if myid==0:
         store_parameters(**kwargs)
-    barrier
+   # barrier
     
     print 'memory usage before del domain1',mem_usage()
     
-def export_model(**kwargs):    
-    #store_parameters(**kwargs)
-    
-#    print 'memory usage before del domain',mem_usage()
-    #del domain
-    print 'memory usage after del domain',mem_usage()
-    
-    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
-    print'swwfile',swwfile
-    
-    export_grid(swwfile, extra_name_out = 'town',
-            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
-            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
-            timestep = None,
-            reduction = max,
-            cellsize = kwargs['export_cellsize'],
-            NODATA_value = -9999,
-            easting_min = project.eastingmin,
-            easting_max = project.eastingmax,
-            northing_min = project.northingmin,
-            northing_max = project.northingmax,
-            verbose = False,
-            origin = None,
-            datum = 'GDA94',
-            format = 'asc')
-            
-    inundation_damage(swwfile+'.sww', project.buildings_filename, 
-                      project.buildings_filename_out)
     
 #-------------------------------------------------------------
@@ -269,10 +262,8 @@
     kwargs['starttime']=project.starttime
     kwargs['yieldstep']=project.yieldstep
-    kwargs['midtime']=project.midtime
     kwargs['finaltime']=project.finaltime
     
     kwargs['output_dir']=project.output_run_time_dir
-    kwargs['bathy_file']=project.combined_dir_name
-#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
+    kwargs['elevation_file']=project.combined_dir_name+'.txt'
     kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
     kwargs['file_name']=project.home+'detail.csv'
@@ -291,5 +282,3 @@
     run_model(**kwargs)
       
-    if myid==0:
-        export_model(**kwargs)
-    barrier
+    #barrier
