source: production/gippsland_2005/get_timeseries.py @ 3330

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

"""

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 *


    
#swwfile = project.newoutputname + '.sww'
swwfile = project.outputname

# Causes a memory error
#swwfile = 'O:/2/cit/inundation/Gippsland Lakes/120106/lakes_100_628759.sww'

swwfile = 'O:/2/cit/inundation/Gippsland Lakes/lakes_100_314920.sww'


#Time interval to plot
tmin = 13000
tmax = 21000

def get_gauges_from_file(filename):
    fid = open(filename)
    lines = fid.readlines()
    fid.close()
    
    gauges = []
    gaugelocation = []
    for line in lines[1:]:
        fields = line.split(',')
        # my gauge file set up as locationname, easting, northing
        location = fields[0]
        easting = float(fields[1])
        northing = float(fields[2])
        #z, easting, northing  = redfearn(lat, lon)
        gauges.append([easting, northing])
        gaugelocation.append(location)
        
    #Return gauges and raw data for subsequent storage
    #return gauges, linesfs
    return gauges, lines, gaugelocation

gauges, lines, locations = get_gauges_from_file(project.gauge_filename)

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


               
from math import sqrt
N = len(gauges)
for k, g in enumerate(gauges):
    if k%((N+10)/10)==0: # diagnostics - print 10 lines
        print 'Doing row %d of %d' %(k, N)

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

    max_depth = 0
    for i, t in enumerate(f.get_time()): # T is a list of times
        #if tmin < t < tmax:
        w = f(t, point_id = k)[0]
        z = f(t, point_id = k)[1]
        #myloc = locations[k]
        
        model_time.append(t)        
        stages.append(w)
        elevations.append(z)  #Should be constant over time

        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)

    plot(model_time, stages, '-b')      
    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')

    

show()