Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/Arc_asc2raster_GDA94z50.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/Arc_asc2raster_GDA94z50.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/Arc_asc2raster_GDA94z50.py	(revision 7678)
@@ -23,10 +23,12 @@
 
 scenario_dir="\\\\nas2\\gemd\\georisk_models\\inundation\\data\\western_australia\\bunbury_storm_surge_scenario_2009\\"
+#scenario_dir = "C:\\WorkSpace\\data\\bunbury_storm_surge_scenario_2009\\"
 output_dir = "anuga\\outputs\\"
 
-time_dir1 = '20100103_102505_run_storm_surge_final_0_alby_coarse_lfountai'
+time_dir1 = '20100311_181323_run_storm_surge_final_0_alby_coarse_fountl' #No storm gate
+time_dir2 = '20100312_132550_run_storm_surge_final_0_alby_coarse_fountl' #Closed storm gate
 
 
-time_dirs = [time_dir1] ##, time_dir2, time_dir3, time_dir4, time_dir5,
+time_dirs = [time_dir1, time_dir2] ##, time_dir2, time_dir3, time_dir4, time_dir5,
              ##time_dir6, time_dir7, time_dir8, time_dir9, time_dir10]
 
@@ -36,15 +38,11 @@
     folder =  scenario_dir + output_dir + time_dir + '\\'
     raster_gbd = folder + 'raster.gdb'
-##    elevation_gdb = folder + 'elevation.gdb'
-    #contour = raster_gbd + '\\contour_dep'
-    land = scenario_dir + "ArcGIS\\data.gdb\\initial_conditions_10m"
-##    ocean = scenario_dir + "map_work\\Perth.gdb\\Outlines\\initial_condition_ocean"
+    land = scenario_dir + "ArcGIS\\data.gdb\\initial_conditions_10m_clip"
+    ocean = scenario_dir + "ArcGIS\\data.gdb\\initial_conditions_ocean"
     
     print 'Process: Create File GDB'
-##    gp.createFileGDB_management(folder, "elevation")
-##    gp.CreateFileGDB_management(folder, "raster")
+    # gp.CreateFileGDB_management(folder, "raster")
 
     gp.Workspace = raster_gbd
-##    gp.Workspace = elevation_gdb
 
     print gp.Workspace
@@ -52,11 +50,14 @@
     #replication dictionary
     replicate = (('gold_coast', ''), ('_time_29220_0', 'b'), ('_time_58440_0', 'c'),
-                 ('_time_28860_0', 'b'), ('_time_57720_0', 'c'),
-                 ('_', ''), ('max','M_'), ('depth','_dep_'),
+                 ('_time_28860_0', 'b'), ('_time_57720_0', 'c'), ('storm_gate', 'gate'),
+                 ('_', ''), ('max','M_'), ('depth','_dep_'), ('albycoarse', ''),
                  ('speed', '_spe_'), ('elevation', '_ele_'), ('stage','_stage'))
 
     generate_filename = []
-##    input_ascii = glob.glob(folder + '*elevation*.asc')
-    input_ascii = glob.glob(folder + '*depth*.asc')
+    # input_ascii = glob.glob(folder + '*elevation*.asc')
+    # input_ascii = glob.glob(folder + '*speed*.asc')
+    # input_ascii = glob.glob(folder + '*depth*.asc')
+    input_ascii = glob.glob(folder + '*stage*.asc')
+    
     print time_dir
 
@@ -86,5 +87,6 @@
     
         print 'Process: Extract by Mask'
-        gp.ExtractByMask_sa(output_DEM, land, output_extract)
+        # gp.ExtractByMask_sa(output_DEM, land, output_extract)
+        # gp.ExtractByMask_sa(output_DEM, ocean, output_extract)
 
        
Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_elevation.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_elevation.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_elevation.py	(revision 7678)
@@ -94,5 +94,7 @@
 G9 = geospatial_data[keylist[8]]
 G10 = geospatial_data[keylist[9]]
-#G11 = geospatial_data[keylist[10]]
+G11 = geospatial_data[keylist[10]]
+G12 = geospatial_data[keylist[11]]
+G13 = geospatial_data[keylist[12]]
 
 
@@ -120,5 +122,5 @@
 #~ print 'G10', G10.attributes
 
-G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 #+ G11
+G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 + G11 + G12 + G13
 
 # G = None
@@ -133,5 +135,5 @@
 
 print 'Export combined DEM file'
-G.export_points_file(project.combined_elevation + '.pts')
+G.export_points_file(project.combined_elevation + '_GEMS.pts')
 
 print 'Do txt version too'
@@ -141,5 +143,5 @@
 	#~ g.export_points_file(join(project.topographies_folder, (str(keylist[i]) +'_export.txt')))
 
-G.export_points_file(project.combined_elevation + '.txt')
+G.export_points_file(project.combined_elevation + '_GEMS.txt')
 
 print 'Export individual text files'
@@ -149,11 +151,13 @@
 G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt'))
 G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt'))
-G3.export_points_file(join(project.topographies_folder, keylist[2] +'_export.txt'))
-G4.export_points_file(join(project.topographies_folder, keylist[3] +'_export.txt'))
-G5.export_points_file(join(project.topographies_folder, keylist[4] +'_export.txt'))
-G6.export_points_file(join(project.topographies_folder, keylist[5] +'_export.txt'))
-G7.export_points_file(join(project.topographies_folder, keylist[6] +'_export.txt'))
-G8.export_points_file(join(project.topographies_folder, keylist[7] +'_export.txt'))
-G9.export_points_file(join(project.topographies_folder, keylist[8] +'_export.txt'))
-G10.export_points_file(join(project.topographies_folder, keylist[9] +'_export.txt'))
-#G11.export_points_file(join(project.topographies_folder, keylist[10] +'_export.txt'))
+G3.export_points_file(join(project.topographies_folder, keylist[2][:-4] +'_export.txt'))
+G4.export_points_file(join(project.topographies_folder, keylist[3][:-4] +'_export.txt'))
+G5.export_points_file(join(project.topographies_folder, keylist[4][:-4] +'_export.txt'))
+G6.export_points_file(join(project.topographies_folder, keylist[5][:-4] +'_export.txt'))
+G7.export_points_file(join(project.topographies_folder, keylist[6][:-4] +'_export.txt'))
+G8.export_points_file(join(project.topographies_folder, keylist[7][:-4] +'_export.txt'))
+G9.export_points_file(join(project.topographies_folder, keylist[8][:-4] +'_export.txt'))
+G10.export_points_file(join(project.topographies_folder, keylist[9][:-4] +'_export.txt'))
+G11.export_points_file(join(project.topographies_folder, keylist[10][:-4] +'_export.txt'))
+G12.export_points_file(join(project.topographies_folder, keylist[11][:-4] +'_export.txt'))
+G13.export_points_file(join(project.topographies_folder, keylist[12][:-4] +'_export.txt'))
Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py	(revision 7678)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py	(revision 7678)
@@ -0,0 +1,145 @@
+# ---------------------------------------------------------------------------
+# This script reads in both the master STS (output from process_gems_grids.py) 
+# and the gems boundary order file (defined as gems_order in project.py)
+# to generate a boundary STS file which contains data just at the boundary 
+# points for all timesteps, which is input into run_model.py in ANUGA.
+# Written by Nariman Habili and Leharne Fountain, 2010
+# ---------------------------------------------------------------------------
+
+import os
+import numpy
+from Scientific.IO.NetCDF import NetCDFFile
+from anuga.utilities.polygon import read_polygon
+import sys
+# sys.path.append('W:\\sandpits\\lfountain\\anuga_work\\production\\mandurah_storm_surge_2009\\gems_comparison')
+sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
+import project_GEMS_grid
+from Scientific.IO.NetCDF import NetCDFFile
+from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference
+
+#-----------------
+# Directory setup
+#-------------------
+state = 'western_australia'
+scenario_folder = 'bunbury_storm_surge_scenario_2009'
+event = 'case_1' # Baseline TC Alby event
+
+ENV_INUNDATIONHOME = 'INUNDATIONHOME'
+home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
+# GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS')
+anuga_folder = os.path.join(home, state, scenario_folder, 'anuga')
+boundaries_folder = os.path.join(anuga_folder, 'boundaries')
+event_folder = os.path.join(boundaries_folder, event)
+
+#--------------------------------------------------------
+# Define boundary points to extract from master STS file
+#--------------------------------------------------------
+gems_boundary = numpy.array(read_polygon(project_GEMS_grid.gems_order))
+number_of_points = gems_boundary.shape[0]
+
+#---------------------------------------------------------
+# Read master sts file generated from GEMS ascii grids
+#---------------------------------------------------------
+print "Reading master NetCDF file"
+infile = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'r')
+
+grid_size = infile.grid_size
+ncols = infile.ncols
+nrows = infile.nrows
+no_data = infile.no_data
+refzone = infile.zone
+x_origin = infile.xllcorner
+y_origin = infile.yllcorner
+x = infile.variables['x'][:]
+y = infile.variables['y'][:]
+time = infile.variables['time'][:]
+elevation = infile.variables['elevation'][:]
+stage_array = infile.variables['stage'][:]
+xmomentum_array = infile.variables['xmomentum'][:]
+ymomentum_array = infile.variables['ymomentum'][:]
+number_of_timesteps = stage_array.shape[0]
+
+infile.close()
+
+#--------------------------------------------
+# Creating index from boundary points
+#--------------------------------------------
+print "Creating index from boundary points" # finding location of boundary points in master file
+
+origin = numpy.array([x_origin[0], y_origin[0]])
+index_array_2D = ((gems_boundary - origin)/grid_size).astype(int)
+assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols
+assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows
+
+index_array_1D = numpy.empty(index_array_2D.shape[0], int)
+index_array_1D = index_array_2D[:,1]*ncols + index_array_2D[:,0]
+
+#--------------------------------------------------
+# Extract values from master sts at index points
+#--------------------------------------------------
+print "Extract values from master sts at index points"
+
+x = x.take(index_array_1D)
+y = y.take(index_array_1D)
+elevation = elevation.take(index_array_1D)
+stage = numpy.empty([number_of_timesteps, number_of_points], float)
+xmomentum = numpy.empty([number_of_timesteps, number_of_points], float)
+ymomentum = numpy.empty([number_of_timesteps, number_of_points], float)
+
+for i, k in enumerate(index_array_1D):
+    stage[:,i] = stage_array[:,k] 
+    xmomentum[:,i] = xmomentum_array[:,k] 
+    ymomentum[:,i] = ymomentum_array[:,k] 
+
+#-----------------------------
+# Create the boundary STS file
+#-----------------------------
+print "Creating the boundary STS NetCDF file"
+
+fid = NetCDFFile(os.path.join(event_folder, event + '_boundary.sts'), 'w')
+fid.institution = 'Geoscience Australia'
+fid.description = "description"
+fid.starttime = 0.0
+fid.ncols = ncols
+fid.nrows = nrows
+fid.grid_size = grid_size
+fid.no_data = no_data
+fid.createDimension('number_of_points', number_of_points)
+fid.createDimension('number_of_timesteps', number_of_timesteps)
+fid.createDimension('numbers_in_range', 2)
+
+fid.createVariable('x', 'd', ('number_of_points',))
+fid.createVariable('y', 'd', ('number_of_points',))
+fid.createVariable('elevation', 'd', ('number_of_points',))
+fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
+fid.createVariable('time', 'i', ('number_of_timesteps',))
+fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
+fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
+fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
+fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
+fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
+fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
+
+fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['elevation'][:] = elevation
+fid.variables['time'][:] = time
+
+s = fid.variables['stage']
+xm = fid.variables['xmomentum']
+ym = fid.variables['ymomentum']
+
+for i in xrange(number_of_timesteps):
+    s[i] = stage[i,:]
+    xm[i] = xmomentum[i,:]
+    ym[i] = ymomentum[i,:]
+ 
+origin = Geo_reference(refzone, x_origin, y_origin)
+geo_ref = write_NetCDF_georeference(origin, fid)
+
+fid.variables['x'][:] = x 
+fid.variables['y'][:] = y 
+
+fid.close()
Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/export_results_max.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/export_results_max.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/export_results_max.py	(revision 7678)
@@ -21,7 +21,8 @@
 
 directory = project.output_folder
-time_dir1 = '20100103_102505_run_storm_surge_final_0_alby_coarse_lfountai'
+time_dir1 = '20100311_181323_run_storm_surge_final_0_alby_coarse_fountl' #No storm gate
+#time_dir2 = '20100312_132550_run_storm_surge_final_0_alby_coarse_fountl' #Closed storm gate
 
-time_dirs = [time_dir1]  
+time_dirs = [time_dir1] #, time_dir2]  
 
 # sww filename extensions ie. hobart_time_17640_0.sww, input into list 17640
@@ -31,10 +32,10 @@
 #Modify the cellsize value to set the size of the raster you require
 #Take into account mesh size when aplying this paramater
-cellsize = 20 #250
+cellsize = 0.5 #10
 
 #Now set the timestep at which you want the raster generated.
 #Either set the actual timestep required or use 'None' to indicate that
 #you want the maximum values in the raster over all timesteps
-timestep = 0
+timestep = None
 
 # Set the special areas of interest.  If none, do: area=['All']
@@ -42,8 +43,8 @@
 
 #area = ['South', 'NW', 'Hobart'] 
-area = ['All']      
+area =  ['Storm_gate'] #['All'] #    
 
 # one or more key strings from var_equations below
-var = ['elevation']
+var = ['depth', 'stage', 'speed']
 
 ######
@@ -100,5 +101,5 @@
                 sww2dem(name, basename_out = outname,
                             quantity = quantityname,
-                            timestep = timestep,
+                            reduction = timestep,
                             cellsize = cellsize,      
                             easting_min = easting_min,
@@ -106,5 +107,4 @@
                             northing_min = northing_min,
                             northing_max = northing_max,        
-                            reduction = max, 
                             verbose = True,
                             format = 'asc')
Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py	(revision 7678)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py	(revision 7678)
@@ -0,0 +1,254 @@
+# ---------------------------------------------------------------------------
+# This python script reads in GEMSURGE outputs and writes them into a master 
+# STS file that contains all the data across the entire model domain at all 
+# timesteps. It requires a vertically flipped elevation ASCII grid which is 
+# merged with the GEMSURGE data to calculate the required quantities of 
+# xmomentum and ymomentum.
+# Written by Nariman Habili and Leharne Fountain, 2010
+# ---------------------------------------------------------------------------
+
+import glob
+import os
+import gzip
+import pdb 
+import numpy
+import math
+from Scientific.IO.NetCDF import NetCDFFile
+from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference
+
+#-----------------
+# Directory setup
+#-------------------
+state = 'western_australia'
+scenario_folder = 'bunbury_storm_surge_scenario_2009'
+event = 'case_1' # Baseline TC Alby event
+
+ENV_INUNDATIONHOME = 'INUNDATIONHOME'
+home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
+# GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS')
+anuga_folder = os.path.join(home, state, scenario_folder, 'anuga')
+boundaries_folder = os.path.join(anuga_folder, 'boundaries')
+event_folder = os.path.join(boundaries_folder, event)
+
+# input_dir = os.path.join(GEMS_folder, event)
+
+#-----------------------
+# Input files from GEMS
+#-----------------------
+elevation_file_name = os.path.join(event_folder, ***'topog_flip.asc'***) # Name of the vertically flipped elevation grid
+grid_file_names = glob.glob(os.path.join(event_folder, '*.gz'))
+grid_file_names.sort()
+number_of_timesteps = len(grid_file_names)
+grid_size = 20 #0.0005 # cellsize from the GEMSURGE output grids
+start_time = 0  # all times in seconds
+end_time = 142560 #88200 # all times in seconds
+timestep = 720 #1800 # all times in seconds
+refzone = 50 # UTM zone
+x_origin = ***#115.550 #364221.03   # xllcorner from GEMSURGE output grids
+y_origin = ***#-32.790 #6371062.71  # yllcorner from GEMSURGE output grids
+event_sts = os.path.join(event_folder, event)
+
+assert number_of_timesteps * timestep == end_time - start_time
+
+elevation = []
+stage = []
+speed = []
+theta = []
+
+elevation_file = open(elevation_file_name, 'rb')
+lines = elevation_file.readlines()
+elevation_file.close()
+
+# Strip off the ASCII header and also read ncols, nrows, 
+# x_origin, y_origin, grid_size, no_data value, and read in elevation grid:
+for i, L in enumerate(lines):
+    if i == 0:
+        ncols = int(L.strip().split()[1])
+    if i == 1:
+        nrows = int(L.strip().split()[1])
+    # if i == 2:
+        # x_origin = int(float(L.strip().split()[1])) # this assumes an integer xllcorner
+    # if i == 3:
+        # y_origin = int(float(L.strip().split()[1])) # this assumes an integer yllcorner
+    # if i == 4:
+        # grid_size = int(float(L.strip().split()[1]))
+    # if i == 5:
+        # no_data = int(float(L.strip().split()[1]))
+    if i > 5:
+        elevation+= L.strip().split()
+
+print 'Number or columns: ', ncols
+print 'Number of rows: ', nrows
+print 'Grid size: ', grid_size
+print 'X origin: ', x_origin
+print 'Y origin: ', y_origin
+  
+for f in grid_file_names:
+    print 'file: ', f
+    
+    gz_file = gzip.open(f, 'rb')
+    lines = gz_file.readlines()
+    gz_file.close()
+
+    for i, L in enumerate(lines):
+        
+        if i > 7 and i < (8 + nrows):        # Number 7 refers to the number of rows of header info that we skip - CHECK format
+            stage += [L.strip().split()]
+            
+        if i > (9 + nrows) and i < (10 + 2*nrows):
+            speed += [L.strip().split()]
+            
+        if i > (11 + 2*nrows):
+            theta += [L.strip().split()]
+
+#------------------------------------------------------
+# Create arrays of elevation, stage, speed and theta
+#-------------------------------------------------------
+print 'creating numpy arrays: elevation, stage, speed, theta'
+            
+elevation = numpy.array(elevation).astype('d')
+
+stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
+speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
+theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
+
+for i in range(number_of_timesteps):
+    ii = nrows*i
+    stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d'))
+    speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d'))
+    theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d'))
+
+stage = stage_array.ravel()
+speed = speed_array.ravel()
+theta = theta_array.ravel()
+    
+print 'Size of elevation array: ', elevation.size
+print 'Size of stage array: ', stage.size
+print 'Size of speed array: ', speed.size
+print 'Size of theta array: ', theta.size
+
+assert stage.size == speed.size == theta.size == ncols * nrows * number_of_timesteps
+assert stage.size == number_of_timesteps * elevation.size
+
+depth = numpy.empty(stage.size, dtype='d')
+number_of_points = ncols * nrows
+
+for i in xrange(number_of_timesteps):
+    j = number_of_points * i
+    k = j + number_of_points
+    depth[j:k] = stage[j:k] - elevation # Check GEMS elevations below sea level are negative
+
+momentum = depth * speed
+
+xmomentum = numpy.empty(momentum.size, dtype='d')
+ymomentum = numpy.empty(momentum.size, dtype='d')
+
+print 'Calculating xmomentum and ymomentum'
+
+for i, t in enumerate(theta):
+
+    if t > 0.0 and t < 90.0:
+        mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
+        my = momentum[i]*math.cos(math.radians(t))
+        assert mx > 0 and my > 0
+    if t > 90.0 and t < 180.0:
+        mx = momentum[i]*math.sin(math.radians(180 - t))
+        my = momentum[i]*math.cos(math.radians(180 - t))
+        assert mx > 0 and my < 0
+    if t > 180.0 and t < 270.0:
+        mx = momentum[i]*math.cos(math.radians(270 - t))
+        my = momentum[i]*math.sin(math.radians(270 - t))
+        assert mx < 0 and my < 0
+    if t > 270.0 and t < 360.0:
+        mx = momentum[i]*math.sin(math.radians(360 - t))
+        my = momentum[i]*math.cos(math.radians(360 - t))
+        assert mx < 0 and my > 0
+    if t == 0.0 or t == 360.0:
+        mx = 0
+        my = momentum[i]
+        assert my > 0
+    if t == 90.0:
+        mx = momentum[i]
+        my = 0
+        assert mx > 0
+    if t == 180.0:
+        mx = 0
+        my = momentum[i]
+        assert my < 0
+    if t == 270.0:
+        mx = momentum[i]
+        my = 0
+        assert mx < 0
+        
+    xmomentum[i] = mx
+    ymomentum[i] = my
+
+    assert mx[i]^2 + my[i]^2 == momentum[i]
+    
+x_origin_int = int(10000*x_origin)
+y_origin_int = int(10000*y_origin)
+grid_size_int = int(10000*grid_size)
+
+x = numpy.tile(numpy.arange(x_origin_int, (x_origin_int + ncols * grid_size_int), grid_size_int)/10000.0, nrows)
+y = numpy.repeat(numpy.arange(y_origin_int, (y_origin_int + nrows * grid_size_int), grid_size_int)/10000.0, ncols)
+
+assert x.size == y.size == number_of_points
+
+time = numpy.arange(start_time, end_time, timestep, dtype='i')
+
+assert time.size == number_of_timesteps
+assert momentum.size == stage.size
+
+#-----------------------------
+# Create the STS file
+#-----------------------------
+print "Creating the STS NetCDF file"
+
+fid = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'w')
+fid.institution = 'Geoscience Australia'
+fid.description = "description"
+fid.starttime = 0.0
+fid.ncols = ncols
+fid.nrows = nrows
+fid.grid_size = grid_size
+fid.no_data = no_data
+fid.createDimension('number_of_points', number_of_points)
+fid.createDimension('number_of_timesteps', number_of_timesteps)
+fid.createDimension('numbers_in_range', 2)
+
+fid.createVariable('x', 'd', ('number_of_points',))
+fid.createVariable('y', 'd', ('number_of_points',))
+fid.createVariable('elevation', 'd', ('number_of_points',))
+fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
+fid.createVariable('time', 'i', ('number_of_timesteps',))
+fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
+fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
+fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
+fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
+fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
+fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
+
+fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
+fid.variables['elevation'][:] = elevation
+fid.variables['time'][:] = time
+
+s = fid.variables['stage']
+xm = fid.variables['xmomentum']
+ym = fid.variables['ymomentum']
+
+for i in xrange(number_of_timesteps):
+    ii = i*number_of_points
+    s[i] = stage[ii : ii + number_of_points]
+    xm[i] = xmomentum[ii : ii + number_of_points]
+    ym[i] = ymomentum[ii : ii + number_of_points]
+	
+origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings - it works for long and lat
+geo_ref = write_NetCDF_georeference(origin, fid)
+
+fid.variables['x'][:] = x - geo_ref.get_xllcorner()
+fid.variables['y'][:] = y - geo_ref.get_yllcorner()
+
+fid.close()
Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/project.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/project.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/project.py	(revision 7678)
@@ -38,5 +38,6 @@
 friction = 0.01           # manning's friction coefficient
 starttime = 86400          # start time for simulation - -equivalent to 0000h 4 April 1978
-finaltime = 172800         # final time for simulation - 24 hours for TC Alby
+duration = 86400
+#finaltime = 172800         # final time for simulation - 24 hours for TC Alby
 
 setup = 'storm_surge_final'      # This can be one of four values
@@ -83,18 +84,22 @@
 # Format for point is x,y,elevation (with header)
 point_filenames = ['DPI.txt',							# Bathymetry data from DPI
-                   'Busselton_Chart_Clip_ss.txt',		# Clipped from Busselton_Chart  - see Busselton Tsunami Scenario 2009
-                   'Busselton_NavyFinal_Clip_ss.txt',	# Clipped from Busselton_NavyFinal - see Busselton Tsunami Scenario 2009
-                   'DPI5U1A02_01a_edited_v2.txt',          # Bathymetric LiDAR from DPI - split into manageable pieces and edited so 
-                   'DPI5U1A02_01b_edited_v2.txt',	        # depths below 0 m are negative, and all soundings on land (ie positive) 
-                   'DPI5U1A02_01c_edited_v2.txt',          # are removed as these are not corrected to "bare earth".
+                   'Busselton_Chart.txt',		        # see Busselton Tsunami Scenario 2009
+                   'Busselton_NavyFinal.txt',	        # see Busselton Tsunami Scenario 2009
+                   'Busselton_250m.txt',                # see Busselton Tsunami Scenario 2009
+                   'DPI5U1A02_01a_edited_v2.txt',       # Bathymetric LiDAR from DPI - split into manageable pieces and edited so 
+                   'DPI5U1A02_01b_edited_v2.txt',	    # depths below 0 m are negative, and all soundings on land (ie positive) 
+                   'DPI5U1A02_01c_edited_v2.txt',       # are removed as these are not corrected to "bare earth".
                    'DPI5U1A02_01d_edited_v2.txt',
-                   'DPI5U1A02_01e_edited_v2.txt']
-#                   'Leschenault_TIN.txt']               # TIN created over the Leschenault Estuary and Inlet]					
+                   'DPI5U1A02_01e_edited_v2.txt',
+                   'Bunbury_2010_shoalest.csv',         # Multibeam infill for areas not captured in the Bathy LiDAR
+                   # 'bu0506hy_cut.csv',                  # Hydro survey of Leschenault Inlet and Estuary, extra colums cut
+                   'Leschenault_TIN.txt']               # TIN created over the Leschenault Estuary and Inlet]					
 
 # BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
 # Used in build_elevation.py & run_model.py
 # Format for points easting,northing (no header)
-bounding_polygon_filename = 'bounding_polygon_ss.csv'
+bounding_polygon_filename = 'bounding_polygon_GEMS.csv' #'bounding_polygon_ss.csv'
 bounding_polygon_maxarea = 50000
+
 
 # INTERIOR REGIONS -  for designing the mesh
@@ -102,7 +107,8 @@
 # Format for points easting,northing (no header)                    
 interior_regions_data = [['intermediate.csv', 2500],
-                         ['area_of_interest.csv', 100],
-                         ['storm_gate_area.csv', 1],
-                         ['stormgates.csv', 1]]
+                         #['aoi_small.csv',200],
+                         ['area_of_interest.csv', 200],
+                         ['storm_gate_area.csv', 20]]
+                         # ['stormgates.csv', 1]]
 
 # LAND - used to set the initial stage/water to be offcoast only
@@ -110,6 +116,6 @@
 land_initial_conditions_filename = [['initial_conditions.csv', 0]]
 
-# GEMS order filename
-# Format is index,northing, easting, elevation (without header)
+# GEMS order filename - should be in same direction as landward boundary points ie clockwise or anti-clockwise
+# Format is northing, easting (without header)
 gems_order_filename = 'gems_boundary_order_thinned.csv'
 
@@ -117,5 +123,5 @@
 # Format is as for a building file to be read by csv2building_polygons, 
 # easting, northing, id, floors (with header)
-storm_gate_filename = 'storm_gates.csv'
+storm_gate_filename = 'storm_gates_large.csv' #'storm_gates.csv'
 
 # GAUGES - for creating timeseries at a specific point
@@ -130,5 +136,5 @@
 
 # Landward bounding points
-# Format easting,northing (no header)
+# Format easting,northing (no header) - should be in same direction as GEMS order filename ie clockwise or anti-clockwise
 landward_boundary_filename = 'landward_boundary.csv'
 
@@ -258,3 +264,10 @@
 ##mux_input = join(event_folder, mux_input_filename)
 
-
+# Areas for export of results
+# Used in export_results_max.py
+
+# Storm gate area:
+xminStorm_gate=373515
+xmaxStorm_gate=373582
+yminStorm_gate=6312314
+ymaxStorm_gate=6312358
Index: /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/run_model.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/run_model.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/run_model.py	(revision 7678)
@@ -38,5 +38,5 @@
 from anuga.interface import Field_boundary
 from anuga.interface import create_sts_boundary
-##from anuga.interface import csv2building_polygons
+from anuga.interface import csv2building_polygons
 from file_length import file_length
 
@@ -45,4 +45,5 @@
 ##from anuga.shallow_water.data_manager import urs2sts
 from anuga.utilities.polygon import read_polygon, Polygon_function
+from anuga.utilities.system_tools import get_revision_number
 
 # Application specific imports
@@ -61,4 +62,6 @@
 start_screen_catcher(project.output_run, 0, 1)
 
+print 'SVN revision number: ', get_revision_number()
+
 #-------------------------------------------------------------------------------
 # Create the computational domain based on overall clipping polygon with
@@ -68,12 +71,4 @@
 
 print 'Create computational domain'
-
-### Create the STS file
-##print 'project.mux_data_folder=%s' % project.mux_data_folder
-##if not os.path.exists(project.event_sts + '.sts'):
-##    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
-
-### Read in boundary from ordered sts file
-##gems_boundary = create_sts_boundary(project.event_sts)
 
 #reading the GEMS boundary points
@@ -139,5 +134,5 @@
 # Add storm surge gate
 #--------------------------
-storm_gate_polygon, storm_gate_height = csv2building_polygons(project.storm_gate, floor_height=4)
+storm_gate_polygon, storm_gate_height = csv2building_polygons(project.storm_gate, floor_height=4.4)
 
 gate = []
@@ -179,5 +174,5 @@
 t0 = time.time()
 for t in domain.evolve(yieldstep=project.yieldstep, 
-                       finaltime=project.finaltime,
+                       duration=project.duration,
                        skip_initial_step=False): 
     print domain.timestepping_statistics()
