source: anuga_work/production/karratha_2005/get_cross_section.py @ 4347

"""Read in sww file, read points along specified cross section and time and plot
"""

import project
from anuga.pyvolution.util import file_function
from anuga.pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn


#swwfile = project.outputdir + project.basename + '.sww'
#swwfile = project.outputname + '_0.0tide' + '.sww'
swwfile = project.outputname + '.sww'
swwfile = project.outputname + '_notsunami.sww'
#swwfile = project.outputdir + 'karratha_100m_12m.sww'


time = 0

#E-W cross section 
start = [468500, 7715177]
end = [469500, 7715177]
step = 1

#Generate points (from west to east)
points = []
x = start[0]
y = start[1]
while x <= end[0]:
    points.append( [x,y] )
    x += step


#Read model output
from caching import cache
f = cache(file_function, swwfile,
          {'quantities': ['stage', 'elevation', 'xmomentum', 'ymomentum'],
           'interpolation_points': points,
           'verbose': True},
          compression = False, 
          dependencies = [swwfile],
          verbose = True)

#print f

T = f.get_time()
from math import sqrt
for time in [T[0], T[0] + 10]:

    abscissa = []
    stages = []
    elevations = []
    momenta = []
    
    for k, p in enumerate(points):
        #
        #print time, k, p
        w = f(time, point_id = k)[0]
        z = f(time, point_id = k)[1]
        uh = f(time, point_id = k)[2]
        vh = f(time, point_id = k)[3]

        abscissa.append( points[k][0] )
        stages.append(w)
        elevations.append(z) 
        momenta.append( sqrt(uh*uh + vh*vh) )  #Absolute momentum

    from pylab import *
    
    ion()
    hold(False)

    if elevations[0] < -10:
        plot(abscissa, stages, '-b')
    else:    
        plot(abscissa, stages, '-b',
             abscissa, elevations, '-k')
    name = 'Cross section (time = %.1f)' %(time)
    title(name)

    title('%s (stage)' %name)
    xlabel('time [s]')
    ylabel('elevation [m]')    
    legend(('Stage', 'Bed = %.1f' %elevations[0]),
           shadow=True,
           loc='upper right')
    savefig('Cross_section_%d' %int(time))

    raw_input('Next')


    #Momentum plot
    ion()
    hold(False)
    plot(abscissa, momenta, '-r')
    title(name)

    title('%s (momentum)' %name)
    xlabel('time [s]')
    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
    savefig('Cross_section_%d_momentum' %int(time))

    raw_input('Next')
    

    
show()