source: anuga_work/development/convergence_okushiri_2008/compare_timeseries.py @ 6409

"""Verify that simulation produced by ANUGA compares to published
validation timeseries ch5, ch7 and ch9 as well as the boundary timeseries.

RMS norm is printed and plots are produced.
"""

from Numeric import allclose, argmin, argmax
from Scientific.IO.NetCDF import NetCDFFile

from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.utilities.numerical_tools import\
     ensure_numeric, cov, get_machine_precision

import project_truescale

try:
    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
except:
    plotting = False
else:
    plotting = True

print 'plotting', plotting
#plotting = False

#-------------------------
# Basic data
#-------------------------

finaltime = 450
timestep = 1

gauge_locations = [[0.000, 678.4]] # Boundary gauge
gauge_locations += [[1808.4, 478.4],  [1808.4, 678.4],  [1808.4, 878.4]] #Ch 5-7-9
gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']

validation_data = {}
for key in gauge_names:
    validation_data[key] = []


#-------------------------
# Read validation dataa
#-------------------------

print 'Reading', project_truescale.boundary_filename
fid = NetCDFFile(project_truescale.boundary_filename, 'r')
input_time = fid.variables['time'][:]
validation_data['Boundary'] = fid.variables['stage'][:]

reference_time = []
fid = open(project_truescale.validation_filename)
lines = fid.readlines()
fid.close()

for i, line in enumerate(lines[1:]):
    if i == len(input_time): break
    
    fields = line.split()

    reference_time.append(float(fields[0]))    # Record reference time
    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
        value = float(fields[1:][j])           # Omit time 
        validation_data[key].append(value) 


# Checks
assert reference_time[0] == 0.0
assert reference_time[-1] == finaltime
assert allclose(reference_time, input_time)

for key in gauge_names:
    validation_data[key] = ensure_numeric(validation_data[key])


#--------------------------------------------------
# Read and interpolate model output
#--------------------------------------------------

import sys
if len(sys.argv) > 1:
    sww_filename = sys.argv[1]
else:    
    sww_filename = project_truescale.output_filename
    
f = file_function(sww_filename,
                  quantities='stage',
                  interpolation_points=gauge_locations,
                  use_cache=True,
                  verbose=True)

#--------------------------------------------------
# Check max runup
#--------------------------------------------------

from anuga.shallow_water.data_manager import get_maximum_inundation_elevation
from anuga.shallow_water.data_manager import get_maximum_inundation_location
from anuga.utilities.polygon import is_inside_polygon

q = get_maximum_inundation_elevation(sww_filename, tolerance=.4)
loc = get_maximum_inundation_location(sww_filename, tolerance=.4)

print 'Max runup elevation: ', q
print 'Max runup location:  ', loc



#--------------------------------------------------
# Compare model output to validation data
#--------------------------------------------------

eps = get_machine_precision()
for k, name in enumerate(gauge_names):
    sqsum = 0
    denom = 0
    model = []
    print 
    print 'Validating ' + name
    observed_timeseries = validation_data[name]
    for i, t in enumerate(reference_time):
        model.append(f(t, point_id=k)[0])


    # Covariance measures    
    res = cov(observed_timeseries, model)     
    print 'Covariance = %.18e' %res

    # Difference measures    
    res = sum(abs(observed_timeseries-model))/len(model)     
    print 'Accumulated difference = %.18e' %res


    # Extrema
    res = abs(max(observed_timeseries)-max(model))
    print 'Difference in maxima = %.18e' %res

    
    res = abs(min(observed_timeseries)-min(model))
    print 'Difference in minima = %.18e' %res

    # Locations of extrema
    i0 = argmax(observed_timeseries)
    i1 = argmax(model)
    res = abs(reference_time[i1] - reference_time[i0])
    print 'Timelag between maxima = %.18e' %res
    

    i0 = argmin(observed_timeseries)
    i1 = argmin(model)
    res = abs(reference_time[i1] - reference_time[i0])
    print 'Timelag between minima = %.18e' %res


    # Write timeseries at gauge locations to file
    filename=name+'_truescale_model.txt'
    fid=open(filename,'w')
    fid.write('Time (s), Stage (m)\n')
    for i, t in enumerate(reference_time):
        v=model[i]
        fid.write('%.2f, %.2f\n' %(t, v))

    fid.close()        
    

    if plotting is True:
        ion()
        hold(False)
    
        plot(reference_time, validation_data[name], 'r-',
             reference_time, model, 'k-')
        title('Gauge %s' %name)
        xlabel('time(s)')
        ylabel('stage (m)')    
        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
        savefig(name, dpi = 300)        

        raw_input('Next')