Index: /production/sydney_2006/export_results.py
===================================================================
--- /production/sydney_2006/export_results.py	(revision 2313)
+++ /production/sydney_2006/export_results.py	(revision 2313)
@@ -0,0 +1,40 @@
+import project, os
+
+from pyvolution.data_manager import sww2dem
+from pyvolution.ermapper_grids import read_ermapper_grid
+
+#name = project.newoutputname
+name = project.outputname
+
+#print 'Which variable do you want to export?'
+#which_var = int(raw_input('Stage = 0, Absolute Momentum = 1, Depth = 2'  ))
+which_var = 2
+
+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'  #Absolute momentum
+
+if which_var == 2:  # Depth
+    outname = name + '_depth'
+    quantityname = 'stage-elevation'  #Depth
+    
+sww2dem(name, basename_out = outname,
+            quantity = quantityname,
+            cellsize = 25,       # Trevor would like this at 25
+            # define region for viz purposes
+            easting_min = project.eminviz,
+            easting_max = project.emaxviz,
+            northing_min = project.nminviz,
+            northing_max = project.nmaxviz,        
+            reduction = max, #this is because we want max quantityname
+            verbose = True,
+            format = 'asc')
+
+#Check
+
+#data = read_ermapper_grid(name)
+#print 'Values from %s are in [%f, %f]' %(name, min(data.flat), max(data.flat))
Index: /production/sydney_2006/get_timeseries.py
===================================================================
--- /production/sydney_2006/get_timeseries.py	(revision 2313)
+++ /production/sydney_2006/get_timeseries.py	(revision 2313)
@@ -0,0 +1,132 @@
+"""Read in sww file, interpolate at specified locations and plot time series
+
+"""
+
+import project
+from pyvolution.util import file_function
+from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+from pylab import *
+
+
+    
+#swwfile = project.newoutputname + '.sww'
+swwfile = project.outputname + '.sww'
+
+#Time interval to plot
+tmin = 13000
+tmax = 21000
+
+def get_gauges_from_file(filename):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    
+    gauges = []
+    gaugelocation = []
+    for line in lines[1:]:
+        fields = line.split(',')
+        # my gauge file set up as locationname, easting, northing
+        location = fields[0]
+        easting = float(fields[1])
+        northing = float(fields[2])
+        #z, easting, northing  = redfearn(lat, lon)
+        gauges.append([easting, northing])
+        gaugelocation.append(location)
+        
+    #Return gauges and raw data for subsequent storage
+    #return gauges, linesfs
+    return gauges, lines, gaugelocation
+
+#gauges, buildings = get_gauges_from_file(project.gauge_filename)
+gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
+
+#Read model output
+quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+f = file_function(swwfile,
+                  quantities = quantities,
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+
+               
+from math import sqrt
+N = len(gauges)
+for k, g in enumerate(gauges):
+    if k%((N+10)/10)==0: # diagnostics - print 10 lines
+        print 'Doing row %d of %d' %(k, N)
+
+    model_time = []
+    stages = []
+    elevations = []
+    momenta = []
+
+    max_depth = 0
+    max_momentum = 0    
+    for i, t in enumerate(f.T): # T is a list of times
+        #if tmin < t < tmax:
+        w = f(t, point_id = k)[0]
+        z = f(t, point_id = k)[1]
+        uh = f(t, point_id = k)[2]
+        vh = f(t, point_id = k)[3]
+        #myloc = locations[k]
+
+        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
+        
+        model_time.append(t)        
+        stages.append(w)
+        elevations.append(z)  #Should be constant over time
+        momenta.append(m) 
+
+        if w-z > max_depth:
+            max_depth = w-z
+        if m > max_momentum:            
+            max_momentum = m    
+
+    
+
+    #Plot only those gauges that have been inundated by more than a threshold
+    #if max_depth < 0.2:
+    #    print 'Skipping gauge %d' %k
+    #    continue
+
+    ion()
+    hold(False)
+
+    if elevations[0] < -10:
+        plot(model_time, stages, '-b')
+    else:    
+        plot(model_time, stages, '-b',
+             model_time, elevations, '-k')
+    name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
+    #name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
+    title(name)
+
+    title('%s (stage)' %name)
+    xlabel('time [s]')
+    ylabel('elevation [m]')    
+    legend(('Stage', 'Bed = %.1f' %elevations[0]),
+           shadow=True,
+           loc='upper right')
+    savefig('Gauge_%d_stage' %k) # savefig('Gauge_%s_stage' %myloc)
+
+    raw_input('Next')
+
+
+    #Momentum plot
+    ion()
+    hold(False)
+    plot(model_time, momenta, '-r')
+    title(name)
+
+    title('%s (momentum)' %name)
+    xlabel('time [s]')
+    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
+    savefig('Gauge_%d_momentum' %k)
+    #savefig('Gauge_%s_momentum' %myloc)
+    
+    raw_input('Next')
+    
+
+show()
+
Index: /production/sydney_2006/process_gauges.py
===================================================================
--- /production/sydney_2006/process_gauges.py	(revision 2313)
+++ /production/sydney_2006/process_gauges.py	(revision 2313)
@@ -0,0 +1,95 @@
+"""Read in sww file, interpolate at specified locations and plot time series
+Additionally, store augmented building database with max inundation depths
+"""
+
+import project
+from pyvolution.util import file_function
+from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+from pylab import *
+
+
+plot = False
+    
+swwfile = project.outputname + '.sww'
+
+def get_gauges_from_file(filename):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+
+    gauges = []
+    gaugelocation = []
+    for line in lines[1:]:
+        fields = line.split(',')
+        location = fields[0]
+        easting = float(fields[1])
+        northing = float(fields[2])
+     
+        gauges.append([easting, northing])
+        gaugelocation.append(location)
+
+    #Return gauges and raw data for subsequent storage    
+    return gauges, lines, gaugelocation
+
+gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
+
+#Read model output
+quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+f = file_function(swwfile,
+                  quantities = quantities,
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+
+lines[0] = lines[0].strip() +\
+               ',MAX INUNDATION DEPTH (m)' +\
+               ',MAX MOMENTUM (m^2/s)\n'            
+               
+from math import sqrt
+N = len(gauges)
+for k, g in enumerate(gauges):
+    if k%((N+10)/10)==0:
+        print 'Doing row %d of %d' %(k, N)
+
+    model_time = []
+    stages = []
+    elevations = []
+    momenta = []
+
+    max_depth = 0
+    max_momentum = 0    
+    for i, t in enumerate(f.T):
+        w = f(t, point_id = k)[0]
+        z = f(t, point_id = k)[1]
+        uh = f(t, point_id = k)[2]
+        vh = f(t, point_id = k)[3]        
+
+        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
+        
+        model_time.append(t)        
+        stages.append(w)
+        elevations.append(z)  #Should be constant over time
+        momenta.append(m) 
+
+        if w-z > max_depth:
+            max_depth = w-z
+        if m > max_momentum:            
+            max_momentum = m    
+
+    #Augment building data
+    lines[k+1] = lines[k+1].strip() +\
+                     ',%f' %max_depth +\
+                     ',%f\n' %max_momentum
+    
+    print lines[k]
+    
+#Store new building file with mak depths added
+
+#FN = 'inundation_augmented_' + project.gauge_filename
+FN = project.gauge_outname
+fid = open(FN, 'w')
+for line in lines:
+    fid.write(line)
+fid.close()    
+
