source: production/karratha_2005/process_buildings.py @ 2898

"""Read in sww file, interpolate at specified locations and plot time series
Additionally, store augmented building database with max inundation depths
"""

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


plot = False
    
swwfile = project.outputname + '.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[4])
        lat = float(fields[5])

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

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

gauges, buildings = get_gauges_from_file(project.gauge_filename)

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

#New header in spreadsheet
buildings[0] = buildings[0].strip() +\
               ',MAX INUNDATION DEPTH (m)' +\
               ',MAX MOMENTUM (m^2/s)\n'            
               
from math import sqrt
N = len(gauges)
for k, g in enumerate(gauges):
    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.get_time()):
        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    

    #Augment building data
    buildings[k+1] = buildings[k+1].strip() +\
                     ',%f' %max_depth +\
                     ',%f\n' %max_momentum
    
    #print buildings[k]
    
    if plot is True:
        #Plot only those gauges that have been inundated by more than a threshold
        if max_depth < 0.2:
            print 'Skipping gauge %d' %k
            continue

        #FIXME Reduce to tmin and tmax
        #i0 = model_timefind
    
        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')
    

#Store new building file with max depths added

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


if plot is True:    
    show()