source: development/okushiri_2005/extract_timeseries.py @ 2230

Last change on this file since 2230 was 2229, checked in by steve, 18 years ago

Moved directories into production and development parent directories

File size: 2.7 KB
Line 
1"""Read in tms file, interpolate at specified locations (ch 5,7,9) and compare
2"""
3
4import sys
5from os import sep
6from caching import cache
7
8sys.path.append('..'+sep+'..'+sep)
9
10
11gauges = [[0.000, 1.696]]  #Boundary
12gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
13
14gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
15
16
17finaltime = 22.5
18timestep = 0.05
19
20#Read reference data
21
22#Input wave
23filename = 'Benchmark_2_input.tms'
24print 'Reading', filename
25from Scientific.IO.NetCDF import NetCDFFile
26from Numeric import allclose
27fid = NetCDFFile(filename, 'r')#
28
29input_time = fid.variables['time'][:]
30input_stage = fid.variables['stage'][:]
31
32
33#gauges
34reference_time = []
35ch5 = []
36ch7 = []
37ch9 = []
38filename = 'output_ch5-7-9.txt'
39fid = open(filename)
40lines = fid.readlines()
41fid.close()
42for i, line in enumerate(lines[1:]):
43    if i == len(input_time): break
44
45    fields = line.split()
46
47    reference_time.append(float(fields[0]))
48    ch5.append(float(fields[1])/100)   #cm2m
49    ch7.append(float(fields[2])/100)   #cm2m
50    ch9.append(float(fields[3])/100)   #cm2m
51
52
53from pyvolution.util import file_function, ensure_numeric
54gauge_values = [ensure_numeric(input_stage),
55                ensure_numeric(ch5),
56                ensure_numeric(ch7),
57                ensure_numeric(ch9)] #Reference values
58
59
60
61#Read model output
62#filename = 'output.sww'
63filename = 'lwru2.sww'
64#filename = 'lwru2_.05s.sww'
65
66#f = file_function(filename,
67#                  quantities = 'stage',
68#                  interpolation_points = gauges,
69#                  verbose = True)
70
71f = cache(file_function,filename,
72          {'quantities': 'stage',
73           'interpolation_points': gauges,
74           'verbose': True},
75          #evaluate = True,
76          dependencies = [filename],
77          verbose = True)
78
79
80
81
82#Checks
83#print reference_time
84#print input_time
85assert reference_time[0] == 0.0
86assert reference_time[-1] == finaltime
87assert allclose(reference_time, input_time)
88
89
90
91#Validation
92
93
94for k, name in enumerate(gauge_names):
95    sqsum = 0
96    denom = 0
97    model = []
98    print 'Validating ' + name
99    for i, t in enumerate(reference_time):
100        ref = gauge_values[k][i]
101        val = f(t, point_id = k)[0]
102        model.append(val)
103
104        sqsum += (ref - val)**2
105        denom += ref**2
106
107    print sqsum
108    print sqsum/denom
109
110    from pylab import *
111    ion()
112    hold(False)
113    plot(reference_time, gauge_values[k], 'r.-',
114         reference_time, model, 'k-')
115    title('Gauge %s' %name)
116    xlabel('time(s)')
117    ylabel('stage (m)')   
118    legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
119    savefig(name, dpi = 300)
120
121    raw_input('Next')
122show()
123
124
125
126#from pylab import *
127#plot(time, stage)
128#show()
Note: See TracBrowser for help on using the repository browser.