1 | """Compare selected timeseries from model outputs |
---|
2 | |
---|
3 | Usage: |
---|
4 | python compare_model_timeseries.py swwfile1 swwfile2 |
---|
5 | |
---|
6 | Return 0 if timeseries are 'close enough', otherwise return non zero error code |
---|
7 | |
---|
8 | Typically swwfile1 would be model output from unit test and swwfile2 would |
---|
9 | be a reference model. |
---|
10 | """ |
---|
11 | |
---|
12 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
13 | import numpy as num |
---|
14 | import sys, os |
---|
15 | |
---|
16 | # Model specific data |
---|
17 | gauges = {'g10a': [422233.4, 874380.2], |
---|
18 | 'g10b': [422107.9, 873556.8], |
---|
19 | 'g10c': [421966.8, 872890.3], |
---|
20 | 'g10d': [421606.1, 872106.2], |
---|
21 | 'g11a': [422628.2, 873933.2], |
---|
22 | 'g11b': [422716.2, 873420.6], |
---|
23 | 'g11c': [422689.1, 872859.8], |
---|
24 | 'g11d': [422408.7, 871940.3], |
---|
25 | 'g11e': [421751.7, 871097.8], |
---|
26 | 'north': [422572.4, 873992.6], |
---|
27 | 'south': [422004.4, 872014.4]} |
---|
28 | |
---|
29 | |
---|
30 | #--------------------------------------- |
---|
31 | |
---|
32 | def usage(): |
---|
33 | print 'Usage:' |
---|
34 | print 'python compare_model_timeseries.py swwfile1 swwfile2' |
---|
35 | |
---|
36 | |
---|
37 | |
---|
38 | def get_timeseries(sww_filename, gauges): |
---|
39 | """Generate time series for sww file based on gauges |
---|
40 | """ |
---|
41 | |
---|
42 | gauge_locations = gauges.values() |
---|
43 | gauge_names = gauges.keys() |
---|
44 | |
---|
45 | print sww_filename |
---|
46 | print gauge_locations |
---|
47 | os.system('cp %s tempfile.sww' % sww_filename) # Has end with sww |
---|
48 | f = file_function('tempfile.sww', #_filename, |
---|
49 | quantities='stage', |
---|
50 | interpolation_points=gauge_locations, |
---|
51 | use_cache=True, |
---|
52 | verbose=True) |
---|
53 | |
---|
54 | timevector = f.get_time() |
---|
55 | |
---|
56 | timeseries = {} |
---|
57 | for k, name in enumerate(gauge_names): |
---|
58 | model = timeseries[name] = [] |
---|
59 | |
---|
60 | for t in timevector: |
---|
61 | model.append(f(t, point_id=k)[0]) |
---|
62 | |
---|
63 | |
---|
64 | return timevector, timeseries |
---|
65 | |
---|
66 | |
---|
67 | |
---|
68 | def compare_timeseries(swwfile1, swwfile2): |
---|
69 | |
---|
70 | timevector1, timeseries1 = get_timeseries(swwfile1, gauges) |
---|
71 | timevector2, timeseries2 = get_timeseries(swwfile2, gauges) |
---|
72 | |
---|
73 | msg = 'Time vectors were different in models %s and %s' %\ |
---|
74 | (swwfile1, swwfile2) |
---|
75 | assert num.allclose(timevector1, timevector2), msg |
---|
76 | |
---|
77 | print 'OK' |
---|
78 | |
---|
79 | |
---|
80 | |
---|
81 | if __name__ == '__main__': |
---|
82 | |
---|
83 | if len(sys.argv) != 3: |
---|
84 | usage() |
---|
85 | import sys; sys.exit(-1) |
---|
86 | |
---|
87 | res = compare_timeseries(sys.argv[1], sys.argv[2]) |
---|
88 | try: |
---|
89 | res = compare_timeseries(sys.argv[1], sys.argv[2]) |
---|
90 | except: |
---|
91 | sys.exit(-1) |
---|
92 | else: |
---|
93 | sys.exit(0) |
---|