Changeset 2069

Timestamp:

Nov 24, 2005, 6:12:22 PM (17 years ago)

Author:

ole

Message:

Karratha work

Files:

• inundation/pyvolution/data_manager.py (modified) (1 diff)
• production/karratha_2005/get_timeseries.py (modified) (5 diffs)
• production/karratha_2005/project.py (modified) (1 diff)

inundation/pyvolution/data_manager.py

@@ -362 +362 @@
         from os import stat
 
-        minimum_allowed_depth = 0.001
-        #minimum_allowed_depth = 0.0  #FIXME pass in or read from domain
+        #minimum_allowed_depth = 0.001
+        minimum_allowed_depth = 0.0  #FIXME pass in or read from domain
         from Numeric import choose
 

production/karratha_2005/get_timeseries.py

@@ -1 @@
-"""Read in sww file, interpolate at specified locations and plot
+"""Read in sww file, interpolate at specified locations and plot time series
+Additionally, store augmented building database with max inundation depths
 """
 
@@ -6 +7 @@
 from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
 from pylab import *
-    
 
 
-#swwfile = project.outputdir + project.basename + '.sww'
-#swwfile = project.outputname + '_0.0tide' + '.sww'
+plot = False
+    
 swwfile = project.outputname + '.sww'
-#swwfile = project.outputname + '_notsunami.sww'
-#swwfile = project.outputdir + 'karratha_100m_12m.sww'
 
 #MOST Boundary directly north of Legendre island
@@ -32 +30 @@
     for line in lines[1:]:
         fields = line.split(',')
-        lon = float(fields[5])
-        lat = float(fields[6])
+        lon = float(fields[4])
+        lat = float(fields[5])
 
         z, easting, northing  = redfearn(lat, lon)        
         gauges.append([easting, northing])
 
-    return gauges
+    #Return gauges and raw data for subsequent storage    
+    return gauges, lines
 
-gauges = get_gauges_from_file(project.gauge_filename)
-
-
-#gauges = gauges[-10:]
-#print gauges
-#import sys; sys.exit() 
-
-#From Neil
-#gauges = [[easting, northing - 5],
-#          [easting, northing - 10],
-#          [easting, northing - 50],
-#          [easting, northing - 100],                                        
-#          [470882, 7717952],
-#          [469390, 7715426],
-#          [469214, 7714979],
-#          [468899, 7715177],
-#          [469038, 7715725],
-#          [470285, 7717470]]
-          
+gauges, buildings = get_gauges_from_file(project.gauge_filename)
 
 #Read model output
@@ -68 +49 @@
                   use_cache = True)
 
-print f
-#model_time = f.T
 
-
+buildings[0] = buildings[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 = []
@@ -80 +65 @@
     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]        
 
-    max_depth = 0
-    for i, t in enumerate(f.T):
-        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]        
+        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) 
 
-            model_time.append(t)        
-            stages.append(w)
-            elevations.append(z)  #Should be constant over time
-            momenta.append( sqrt(uh*uh + vh*vh) )  #Absolute momentum
+        if w-z > max_depth:
+            max_depth = w-z
+        if m > max_momentum:            
+            max_momentum = m    
 
-            if w-z > max_depth:
-                max_depth = w-z
+    #Augment building data
+    buildings[k+1] = buildings[k+1].strip() +\
+                     ',%f' %max_depth +\
+                     ',%f\n' %max_momentum
+    
+    #print buildings[k]
+    
+    if plot is True:
+        #Plot only those gauges that have been inundated by more than a threshold
+        if max_depth < 0.2:
+            print 'Skipping gauge %d' %k
+            continue
 
-                
+        #FIXME Reduce to tmin and tmax
+        #i0 = model_timefind
+    
+        ion()
+        hold(False)
 
-    #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])
+        title(name)
 
-    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)
+        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)
 
-    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')
+        raw_input('Next')
 
 
-    #Momentum plot
-    ion()
-    hold(False)
-    plot(model_time, momenta, '-r')
-    title(name)
+        #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)
+        title('%s (momentum)' %name)
+        xlabel('time [s]')
+        ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
+        savefig('Gauge_%d_momentum' %k)
 
-    raw_input('Next')
+        raw_input('Next')
     
 
-    
-show()
+#Store new building file with mak depths added
 
+FN = 'inundation_augmented_' + project.gauge_filename
+fid = open(FN, 'w')
+for line in buildings:
+    fid.write(line)
+fid.close()    
+
+
+if plot is True:    
+    show()
+

production/karratha_2005/project.py

@@ -31 +31 @@
 outputname = outputdir + basename  #Used by post processing
 
-gauge_filename = 'all_bld_ind.csv'
+#gauge_filename = 'all_bld_ind.csv'
+gauge_filename = 'karratha_buildings_23112005.csv'