Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/Arc_asc2raster_GDA94z50.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/Arc_asc2raster_GDA94z50.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/Arc_asc2raster_GDA94z50.py	(revision 7677)
@@ -25,9 +25,8 @@
 output_dir = "anuga\\outputs\\"
 
-time_dir1 = '20091229_100831_run_storm_surge_final_0_alby_waves_lfountai'
+time_dir1 = '20100219_094651_run_storm_surge_final_0_alby_waves_lfountai'
+# time_dir2 = '20100223_153032_run_storm_surge_final_0_alby_waves_GEMS_compare_lfountai'
 
-
-time_dirs = [time_dir1] ##, time_dir2, time_dir3, time_dir4, time_dir5,
-             ##time_dir6, time_dir7, time_dir8, time_dir9, time_dir10]
+time_dirs = [time_dir1] #, time_dir2] 
 
 for time_dir in time_dirs:
@@ -36,20 +35,15 @@
     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
     
     #replication dictionary
-    replicate = (('gold_coast', ''), ('_time_29220_0', 'b'), ('_time_58440_0', 'c'),
+    replicate = (('alby_waves', 'man'), ('_time_29220_0', 'b'), ('_time_58440_0', 'c'),
                  ('_time_28860_0', 'b'), ('_time_57720_0', 'c'),
                  ('_', ''), ('max','M_'), ('depth','_dep_'),
@@ -57,6 +51,6 @@
 
     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 + '*stage.asc')
     print time_dir
 
@@ -86,5 +80,6 @@
     
         print 'Process: Extract by Mask'
-        gp.ExtractByMask_sa(output_DEM, land, output_extract)
+        # gp.ExtractByMask_sa(output_DEM, land, output_extract) # use for depth to get onshore inundation depths
+        # gp.ExtractByMask_sa(output_DEM, ocean, output_extract) # use for stage to get offshore wave heights
 
        
Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/build_elevation.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/build_elevation.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/build_elevation.py	(revision 7677)
@@ -85,12 +85,12 @@
 # G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box)
 
-G1 = geospatial_data['m_peel_aoi'].clip(project.elevation_clip_box)
-G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box)
-G3 = geospatial_data['DpiU1A03_ss.csv']
-G4 = geospatial_data['MA-46893-SNDS_AHD.csv']
-G5 = geospatial_data['MS0205HY_AHD.csv']
-G6 = geospatial_data['MS0404_AHD.csv']
-G7 = geospatial_data['YU0403HY_AHD.csv']
-G8 = geospatial_data['original_data_ss.csv']
+G1 = geospatial_data[keylist[0]].clip(project.elevation_clip_box)
+G2 = geospatial_data[keylist[1]].clip_outside(project.elevation_clip_box)
+G3 = geospatial_data[keylist[2]]
+G4 = geospatial_data[keylist[3]]
+G5 = geospatial_data[keylist[4]]
+G6 = geospatial_data[keylist[5]]
+G7 = geospatial_data[keylist[6]]
+G8 = geospatial_data[keylist[7]]
 
 G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8
@@ -102,4 +102,4 @@
 # Use for comparision in ARC
 G.export_points_file(project.combined_elevation + '.txt')
-G1.export_points_file(join(project.topographies_folder, (str(project.ascii_grid_filenames[0]) +'_export.txt')))
-G2.export_points_file(join(project.topographies_folder, (str(project.ascii_grid_filenames[1]) + '_export.txt')))
+G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt'))
+G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt'))
Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/build_gems_boundary.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/build_gems_boundary.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/build_gems_boundary.py	(revision 7677)
@@ -0,0 +1,157 @@
+# ---------------------------------------------------------------------------
+# 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 = 'mandurah_storm_surge_scenario_2009'
+event = 'Man_1_TCAlby'
+
+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_geo.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"
+
+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('permutation', 'i', ('number_of_thinned_points',))
+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',))
+
+
+#permutation = fid.variables['permutation']
+#elevation = fid.variables['elevation']
+#time = fid.variables['time']
+#stage = fid.variables['stage']
+#xmomentum = fid.variables['xmomentum']
+#ymomentum = fid.variables['ymomentum']
+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
+
+#fid.variables['stage'][:] = stage
+#fid.variables['xmomentum'] = xmomentum
+#fid.variables['ymomentum'] = ymomentum
+
+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[i,:]
+    xm[i] = xmomentum[i,:]
+    ym[i] = ymomentum[i,:]
+ 
+# for i in xrange(number_of_timesteps):#range(1, number_of_points + 1):
+    # _i = i*number_of_points
+    # for j in xrange(number_of_points):
+        # s[j, i] = stage[_i + j]
+        # xm[j, i] = xmomentum[_i + j]
+        # ym[j, i] = ymomentum[_i + j]
+	
+origin = Geo_reference(refzone, x_origin, y_origin)
+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/mandurah_storm_surge_2009/export_results_max.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/export_results_max.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/export_results_max.py	(revision 7677)
@@ -21,7 +21,8 @@
 
 directory = project.output_folder
-time_dir1 = '20091229_100831_run_storm_surge_final_0_alby_waves_lfountai'
+#time_dir1 = '20100219_094651_run_storm_surge_final_0_alby_waves_lfountai'
+time_dir1 = '20100223_153032_run_storm_surge_final_0_alby_waves_GEMS_compare_lfountai'
 
-time_dirs = [time_dir1]  
+time_dirs = [time_dir1] #, time_dir2]  
 
 # sww filename extensions ie. hobart_time_17640_0.sww, input into list 17640
@@ -45,5 +46,5 @@
 
 # one or more key strings from var_equations below
-var = ['stage']
+var = ['elevation', 'speed'] #'stage','depth','speed','elevation']
 
 ######
@@ -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/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py	(revision 7677)
@@ -82,6 +82,6 @@
 print 'project_GEMS.elevation_clip_box', project_GEMS.elevation_clip_box
 
-G1 = geospatial_data['man_25m_prj'].clip(project_GEMS.elevation_clip_box)
-G2 = geospatial_data['man10m_ss'].clip_outside(project_GEMS.elevation_clip_box)
+G1 = geospatial_data[keylist[0]].clip(project.elevation_clip_box)
+G2 = geospatial_data[keylist[1]].clip_outside(project.elevation_clip_box)
 
 G = G1 + G2
@@ -93,4 +93,4 @@
 # Use for comparision in ARC
 G.export_points_file(project_GEMS.combined_elevation + '.txt')
-G1.export_points_file(join(project_GEMS.topographies_folder, (str(project_GEMS.ascii_grid_filenames[0]) +'_export.txt')))
-G2.export_points_file(join(project_GEMS.topographies_folder, (str(project_GEMS.ascii_grid_filenames[1]) + '_export.txt')))
+# G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt'))
+# G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt'))
Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py	(revision 7677)
@@ -180,5 +180,5 @@
 # Output filename for elevation
 # this is a combination of all the data generated in build_elevation.py
-combined_elevation_basename = scenario_name + 'GEMS_combined_elevation'
+combined_elevation_basename = scenario_name + '_GEMS_combined_elevation'
 
 #-------------------------------------------------------------------------------
Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py	(revision 7677)
@@ -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
@@ -60,4 +61,6 @@
                              project_GEMS.__name__+'.py'))
 start_screen_catcher(project_GEMS.output_run, 0, 1)
+
+print 'SVN revision number: ', get_revision_number()
 
 #-------------------------------------------------------------------------------
Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/process_gems_grids.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/process_gems_grids.py	(revision 7677)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/process_gems_grids.py	(revision 7677)
@@ -0,0 +1,243 @@
+# ---------------------------------------------------------------------------
+# 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 = 'mandurah_storm_surge_scenario_2009'
+event = 'Man_1_TCAlby'
+
+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 = 0.0005  # cellsize from the GEMSURGE output grids
+start_time = 0      # all times in seconds
+end_time = 88200    # all times in seconds
+timestep = 1800     # all times in seconds
+refzone = 50        # UTM zone of model
+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)
+
+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, 
+# 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 == 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
+        
+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):   #  7 refers to the number of rows of header info that we skip - always 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 + '_test.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/mandurah_storm_surge_2009/project.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/project.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/project.py	(revision 7677)
@@ -128,5 +128,5 @@
 
 # GEMS order filename
-# Format is index,northing, easting, elevation (without header)
+# Format is index, longitude, latitude, elevation (without header)
 gems_order_filename = 'gems_boundary_order_thinned.csv'
 
Index: /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/run_model.py
===================================================================
--- /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/run_model.py	(revision 7676)
+++ /branches/anuga_1_1/anuga_work/production/mandurah_storm_surge_2009/run_model.py	(revision 7677)
@@ -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
@@ -60,4 +61,6 @@
                              project.__name__+'.py'))
 start_screen_catcher(project.output_run, 0, 1)
+
+print 'SVN revision number: ', get_revision_number()
 
 #-------------------------------------------------------------------------------
