source: anuga_work/development/momentum_sink/scripts/plot_gauges.py @ 6939

"""Read in sww file, interpolate at specified locations and plot time series
Additionally, store augmented building database with max inundation depths
"""
 # FORMALLY CALLED <process_gauges.py>
from os import sep
import project
from anuga.pyvolution.util import file_function
#from anuga.pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
from pylab import *
from compare_sww import gauge_locations

plot = False
    
swwfile = project.outputdir + sep  + 'Buildings_3662.sww'

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

#    gauges = []
#    gaugelocation = []
#    for line in lines[1:]:
#        fields = line.split(',')
#        location = fields[0]
#        easting = float(fields[1])
#        northing = float(fields[2])
     
#        gauges.append([easting, northing])
#        gaugelocation.append(location)

#    #Return gauges and raw data for subsequent storage    
#    return gauges, lines, gaugelocation

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


#Read model output
quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
f = file_function(swwfile,
                  quantities = quantities,
                  interpolation_points = gauge_locations,
                  verbose = True,
                  use_cache = True)     
               
from math import sqrt
N = len(gauge_locations)
for k, g in enumerate(gauge_locations):
    if k%((N+10)/10)==0:
        print 'Doing row %d of %d' %(k, N)

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

    max_depth = 0
    max_momentum = 0    
    for i, t in enumerate(f.T):
        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]        

        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
        
        model_time.append(t)        
        stages.append(w)
        elevations.append(z)  #Should be constant over time
        momenta.append(m) 

        if w-z > max_depth:
            max_depth = w-z
        if m > max_momentum:            
            max_momentum = m    

        plot(model_time, stages, 'b')
    
    print 
    
#Store new building file with mak depths added

#FN = 'inundation_augmented_' + project.gauge_filename
#FN = project.gauge_outname
#fid = open(FN, 'w')
#for line in lines:
#    fid.write(line)
#fid.close()