source: development/momentum_sink/scripts/plot_cross_section.py @ 3322

"""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()