source: development/okushiri_2005/extract_timeseries.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 17 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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 anuga.pyvolution.util import file_function
54from anuga.utilities.numerical_tools import ensure_numeric
55gauge_values = [ensure_numeric(input_stage),
56                ensure_numeric(ch5),
57                ensure_numeric(ch7),
58                ensure_numeric(ch9)] #Reference values
59
60
61
62#Read model output
63#filename = 'output.sww'
64filename = 'lwru2.sww'
65#filename = 'lwru2_.05s.sww'
66
67#f = file_function(filename,
68#                  quantities = 'stage',
69#                  interpolation_points = gauges,
70#                  verbose = True)
71
72f = cache(file_function,filename,
73          {'quantities': 'stage',
74           'interpolation_points': gauges,
75           'verbose': True},
76          #evaluate = True,
77          dependencies = [filename],
78          verbose = True)
79
80
81
82
83#Checks
84#print reference_time
85#print input_time
86assert reference_time[0] == 0.0
87assert reference_time[-1] == finaltime
88assert allclose(reference_time, input_time)
89
90
91
92#Validation
93
94
95for k, name in enumerate(gauge_names):
96    sqsum = 0
97    denom = 0
98    model = []
99    print 'Validating ' + name
100    for i, t in enumerate(reference_time):
101        ref = gauge_values[k][i]
102        val = f(t, point_id = k)[0]
103        model.append(val)
104
105        sqsum += (ref - val)**2
106        denom += ref**2
107
108    print sqsum
109    print sqsum/denom
110
111    from pylab import *
112    ion()
113    hold(False)
114    plot(reference_time, gauge_values[k], 'r.-',
115         reference_time, model, 'k-')
116    title('Gauge %s' %name)
117    xlabel('time(s)')
118    ylabel('stage (m)')   
119    legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
120    savefig(name, dpi = 300)
121
122    raw_input('Next')
123show()
124
125
126
127#from pylab import *
128#plot(time, stage)
129#show()
Note: See TracBrowser for help on using the repository browser.