source: development/stochastic_study/extract_timeseries.py @ 2610

Last change on this file since 2610 was 2610, checked in by ole, 18 years ago

Fixed import

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