1 | """Read in sww file, interpolate at specified locations and plot time series |
---|
2 | |
---|
3 | """ |
---|
4 | |
---|
5 | import project |
---|
6 | from anuga.pyvolution.util import file_function |
---|
7 | from anuga.pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn |
---|
8 | from pylab import * |
---|
9 | |
---|
10 | |
---|
11 | |
---|
12 | swwfile = project.outputname + '.sww' |
---|
13 | |
---|
14 | #MOST Boundary directly north of Legendre island |
---|
15 | lat = project.north |
---|
16 | lon = (project.p3[1] + project.p4[1])/2 |
---|
17 | z, easting, northing = redfearn(lat, lon) |
---|
18 | |
---|
19 | #Time interval to plot |
---|
20 | tmin = 13000 |
---|
21 | tmax = 21000 |
---|
22 | |
---|
23 | def get_gauges_from_file(filename): |
---|
24 | fid = open(filename) |
---|
25 | lines = fid.readlines() |
---|
26 | fid.close() |
---|
27 | |
---|
28 | gauges = [] |
---|
29 | for line in lines[1:]: |
---|
30 | fields = line.split(',') |
---|
31 | lon = float(fields[4]) |
---|
32 | lat = float(fields[5]) |
---|
33 | |
---|
34 | z, easting, northing = redfearn(lat, lon) |
---|
35 | gauges.append([easting, northing]) |
---|
36 | |
---|
37 | #Return gauges and raw data for subsequent storage |
---|
38 | return gauges, lines |
---|
39 | |
---|
40 | gauges, buildings = get_gauges_from_file(project.gauge_filename) |
---|
41 | |
---|
42 | gauges = [[easting, northing-10]] |
---|
43 | |
---|
44 | |
---|
45 | #Read model output |
---|
46 | quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
47 | f = file_function(swwfile, |
---|
48 | quantities = quantities, |
---|
49 | interpolation_points = gauges, |
---|
50 | verbose = True, |
---|
51 | use_cache = True) |
---|
52 | |
---|
53 | |
---|
54 | |
---|
55 | from math import sqrt |
---|
56 | N = len(gauges) |
---|
57 | for k, g in enumerate(gauges): |
---|
58 | if k%((N+10)/10)==0: |
---|
59 | print 'Doing row %d of %d' %(k, N) |
---|
60 | |
---|
61 | model_time = [] |
---|
62 | stages = [] |
---|
63 | elevations = [] |
---|
64 | momenta = [] |
---|
65 | |
---|
66 | max_depth = 0 |
---|
67 | max_momentum = 0 |
---|
68 | for i, t in enumerate(f.get_time()): # T is a list of times |
---|
69 | #if tmin < t < tmax: |
---|
70 | w = f(t, point_id = k)[0] |
---|
71 | z = f(t, point_id = k)[1] |
---|
72 | uh = f(t, point_id = k)[2] |
---|
73 | vh = f(t, point_id = k)[3] |
---|
74 | |
---|
75 | m = sqrt(uh*uh + vh*vh) #Absolute momentum |
---|
76 | |
---|
77 | model_time.append(t) |
---|
78 | stages.append(w) |
---|
79 | elevations.append(z) #Should be constant over time |
---|
80 | momenta.append(m) |
---|
81 | |
---|
82 | if w-z > max_depth: |
---|
83 | max_depth = w-z |
---|
84 | if m > max_momentum: |
---|
85 | max_momentum = m |
---|
86 | |
---|
87 | |
---|
88 | |
---|
89 | #Plot only those gauges that have been inundated by more than a threshold |
---|
90 | #if max_depth < 0.2: |
---|
91 | # print 'Skipping gauge %d' %k |
---|
92 | # continue |
---|
93 | |
---|
94 | ion() |
---|
95 | hold(False) |
---|
96 | |
---|
97 | if elevations[0] < -10: |
---|
98 | plot(model_time, stages, '-b') |
---|
99 | else: |
---|
100 | plot(model_time, stages, '-b', |
---|
101 | model_time, elevations, '-k') |
---|
102 | name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1]) |
---|
103 | title(name) |
---|
104 | |
---|
105 | title('%s (stage)' %name) |
---|
106 | xlabel('time [s]') |
---|
107 | ylabel('elevation [m]') |
---|
108 | legend(('Stage', 'Bed = %.1f' %elevations[0]), |
---|
109 | shadow=True, |
---|
110 | loc='upper right') |
---|
111 | savefig('Gauge_%d_stage' %k) |
---|
112 | |
---|
113 | raw_input('Next') |
---|
114 | |
---|
115 | |
---|
116 | #Momentum plot |
---|
117 | ion() |
---|
118 | hold(False) |
---|
119 | plot(model_time, momenta, '-r') |
---|
120 | title(name) |
---|
121 | |
---|
122 | title('%s (momentum)' %name) |
---|
123 | xlabel('time [s]') |
---|
124 | ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]') |
---|
125 | savefig('Gauge_%d_momentum' %k) |
---|
126 | |
---|
127 | raw_input('Next') |
---|
128 | |
---|
129 | |
---|
130 | show() |
---|
131 | |
---|