source: production/karratha_2005/get_timeseries.py @ 2058

"""Read in sww file, interpolate at specified locations and plot
"""

import project
from pyvolution.util import file_function
from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
from pylab import *
    


#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'

#MOST Boundary directly north of Legendre island
lat = project.north
lon = (project.p3[1] + project.p4[1])/2 
z, easting, northing  = redfearn(lat, lon)

#Time interval to plot
tmin = 13000
tmax = 21000

def get_gauges_from_file(filename):
    fid = open(filename)
    lines = fid.readlines()
    fid.close()

    gauges = []
    for line in lines[1:]:
        fields = line.split(',')
        lon = float(fields[5])
        lat = float(fields[6])

        z, easting, northing  = redfearn(lat, lon)        
        gauges.append([easting, northing])

    return gauges

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]]
          

#Read model output
quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
f = file_function(swwfile,
                  quantities = quantities,
                  interpolation_points = gauges,
                  verbose = True,
                  use_cache = True)

print f
#model_time = f.T


from math import sqrt
for k, g in enumerate(gauges):

    model_time = []
    stages = []
    elevations = []
    momenta = []


    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]        

            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

                

    #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)

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