source: anuga_validation/automated_validation_tests/patong_beach_validation/compare_model_timeseries.py @ 7641

"""Compare selected timeseries from model outputs

Usage:
python compare_model_timeseries.py swwfile1 swwfile2

Return 0 if timeseries are 'close enough', otherwise return non zero error code

Typically swwfile1 would be model output from unit test and swwfile2 would
be a reference model.
"""

import numpy as num
import sys

# Model specific data
gauges = {'g10a': [422233.4, 874380.2],
          'g10b': [422107.9, 873556.8],
          'g10c': [421966.8, 872890.3],
          'g10d': [421606.1, 872106.2],
          'g11a': [422628.2, 873933.2],
          'g11b': [422716.2, 873420.6],
          'g11c': [422689.1, 872859.8],
          'g11d': [422408.7, 871940.3],
          'g11e': [421751.7, 871097.8],
          'north': [422572.4, 873992.6],
          'south': [422004.4, 872014.4]}


#---------------------------------------

def usage():
    print 'Usage:'
    print 'python compare_model_timeseries.py swwfile1 swwfile2'



def get_timeseries(sww_filename, gauges):
    """Generate time series for sww file based on gauges
    """

    gauge_locations = gauges.values()
    gauge_names = gauges.keys()
    
    f = file_function(sww_filename,
                      quantities='stage',
                      interpolation_points=gauge_locations,
                      use_cache=True,
                      verbose=True)

    timevector = f.get_time()

    timeseries = {}
    for k, name in enumerate(gauge_names):
        model = timeseries[name] = []
        
        for t in timevector:
            model.append(f(t, point_id=k)[0])


    return timevector, timeseries



def compare_timeseries(swwfile1, swwfile2):
    
    timevector1, timeseries1 = get_timeseries(swwfile1, gauges)
    timevector2, timeseries2 = get_timeseries(swwfile2, gauges)    

    msg = 'Time vectors were different in models %s and %s' %\
          (swwfile1, swwfile2)
    assert num.allclose(timevector1, timevector2), msg

    


if __name__ == '__main__':

    if len(sys.argv) != 3:
        usage()
        import sys; sys.exit(-1) 
    
    try:
        res = compare_timeseries(sys.argv[1], sys.argv[2])
    except:
        sys.exit(-1)
    else:
        sys.exit(0)