"""Read in sww file, interpolate at specified locations and plot time series """ from os import sep import Numeric import project from pyvolution.util import file_function #from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn #from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn from pylab import * #from compare_sww import gauge_locations #swwfile = project.newoutputname + '.sww' #swwfile = project.outputname swwfile = project.outputdir + sep + 'Buildings_3662.sww' gauge_depth = Numeric.arrayrange(0, 700, 20) gauge_breadth = 100 gauge_locations = [] for GD in gauge_depth: gauge_location = [GD,gauge_breadth] gauge_locations.append(gauge_location) #Read model output quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] f = file_function(swwfile, quantities = quantities, interpolation_points = gauge_locations, verbose = True, use_cache = True) T=[150] from math import sqrt N = len(gauge_locations) for k, g in enumerate(gauge_locations): #if k%((N+10)/10)==0: # diagnostics - print 10 lines # print 'Doing row %d of %d' %(k, N) model_time = [] stages = [] elevations = [] momenta = [] velocity = [] max_depth = 0 max_momentum = 0 max_velocity = 0 for t in T: #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 vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity print vel #dep = w-z #vel = sqrt(uh*uh + vh*vh) / dep #Absolute velocity model_time.append(t) stages.append(w) elevations.append(z) #Should be constant over time momenta.append(m) velocity.append(vel) if w-z > max_depth: max_depth = w-z if m > max_momentum: max_momentum = m if vel > max_velocity: max_velocity = vel print 'max speed', max_velocity print len(stages), "STAGES" #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') # plot(stages, elevations, '-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') #Speed plot #ion() #hold(False) #plot(model_time, velocity, '-r') #title(name) #title('%s (velocity)' %name) #xlabel('time [s]') #ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]') #savefig('Gauge_%d_speed' %k) # savefig('Gauge_%s_speed' %myloc) # raw_input('Next') show()