source: anuga_validation/okushiri_2005/compare_timeseries.py @ 4219

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

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

#plotting = False

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

finaltime = 22.5
timestep = 0.05

gauge_locations = [[0.000, 1.696]] # Boundary gauge
gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #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.boundary_filename
fid = NetCDFFile(project.boundary_filename, 'r')
input_time = fid.variables['time'][:]
validation_data['Boundary'] = fid.variables['stage'][:]

reference_time = []
fid = open(project.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/100) # Convert cm2m


# 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.output_filename
    
f = file_function(sww_filename,
                  quantities='stage',
                  interpolation_points=gauge_locations,
                  use_cache=True,
                  verbose=True)


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




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