Changeset 1921

Timestamp:

Oct 14, 2005, 11:49:24 AM (19 years ago)

Author:

ole

Message:

Get timeseries - work for Karratha

File:

• production/karratha_2005/get_timeseries.py (modified) (4 diffs)

production/karratha_2005/get_timeseries.py

@@ -7 +7 @@
 
 #swwfile = project.outputdir + project.basename + '.sww'
-swwfile = project.outputname + '_0.0tide' + '.sww'
+#swwfile = project.outputname + '_0.0tide' + '.sww'
+swwfile = project.outputname + '.sww'
+#swwfile = project.outputname + '_0.0tide_notsunami.sww'
 #swwfile = project.outputdir + 'karratha_100m_12m.sww'
 
 #MOST Boundary directly north of Legendre island
-y = project.north
-x = (project.p3[1] + project.p4[1])/2   #?????//
-z, easting, northing  = redfearn(y, x)
+lat = project.north
+lon = (project.p3[1] + project.p4[1])/2 
+z, easting, northing  = redfearn(lat, lon)
 
 
 #From Neil
-gauges = [[easting, northing - 10],
+gauges = [[easting, northing - 5],
+          [easting, northing - 10],
+          [easting, northing - 20],
+          [easting, northing - 30],
+          [easting, northing - 50],
+          [easting, northing - 100],                                        
           [470882, 7717952],
           [469390, 7715426],
@@ -25 +32 @@
           [470285, 7717470]]
           
-
-    
-#gauge_names = ['MOST',
-#               'Gauge_1', 'Gauge_2', 'Gauge_3',
-#               'Gauge_4', 'Gauge_5', 'Gauge_6']
- 
 
 #Read model output
@@ -40 +41 @@
 
 f = cache(file_function, swwfile,
-          {'quantities': 'stage',
+          {'quantities': ['stage', 'elevation', 'xmomentum', 'ymomentum'],
            'interpolation_points': gauges,
            'verbose': True},
@@ -51 +52 @@
 
 
+from math import sqrt
+for k, g in enumerate(gauges):
 
-for k, g in enumerate(gauges):
-    model_values = []
+    stages = []
+    elevations = []
+    momenta = []
+    
     for i, t in enumerate(f.T):
-        val = f(t, point_id = k)[0]
-        model_values.append(val)
+        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]        
+        
+        stages.append(w)
+        elevations.append(z)  #Should be constant over time
+        momenta.append( sqrt(uh*uh + vh*vh) )  #Absolute momentum
 
     from pylab import *
+    
     ion()
     hold(False)
-    plot(model_time, model_values, '-b')
-    name = 'Gauge_%d: %s' %(k, str(g))
+
+    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])
     title(name)
-    
-    savefig('Gauge_%d' %k)
+
+    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)
 
     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)
+
+    raw_input('Next')
+    
+
+    
 show()