source: anuga_core/source/anuga/examples/island_timeseries.py @ 3740

Last change on this file since 3740 was 3514, checked in by duncan, 18 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: 5.1 KB
Line 
1"""Read in sww file, interpolate at specified locations and plot time series
2
3"""
4
5from anuga.pyvolution.util import file_function
6from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
7from pylab import *
8from matplotlib.ticker import MultipleLocator, FormatStrFormatter
9
10
11swwfile = 'island.sww' 
12#gauges = [[56, 48], [58, 48], [60, 48], [62, 48], [64, 48]]
13#gauges = [[10, 10], [60, 48], [62, 48], [64, 48]]   
14gauges = [[10, 10], [40, 10], [70, 10]]   
15
16
17store = True
18
19#Read model output
20quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
21f = file_function(swwfile,
22                  quantities = quantities,
23                  interpolation_points = gauges,
24                  verbose = True,
25                  use_cache = True)
26
27
28T = f.get_time()
29               
30from math import sqrt, atan, degrees
31from Numeric import ones
32N = len(gauges)
33for k, g in enumerate(gauges):
34    if k%((N+10)/10)==0: # diagnostics - print 10 lines
35        print 'Doing row %d of %d' %(k, N)
36
37    model_time = []
38    stages = []
39    elevations = []
40    momenta = []
41    velocity = []
42    xmom = []
43    ymom = []
44    bearings = []
45    depths = []
46
47    max_depth = 0
48    max_momentum = 0
49    max_velocity = 0
50
51    due_east = 90.0*ones([len(T)])
52    due_west = 270.0*ones([len(T)])
53    maxT = max(T)
54    tstep = maxT/(len(T)-1)
55
56
57    myloc = 'G%d' %k #locations[k]
58
59    for i, t in enumerate(T): # T is a list of times
60        #if tmin < t < tmax:
61        w = f(t, point_id = k)[0]
62        z = f(t, point_id = k)[1]
63        uh = f(t, point_id = k)[2]
64        vh = f(t, point_id = k)[3]
65        depth = w-z
66
67        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
68        vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
69        angle = degrees(atan(vh/(uh+1.e-15)))
70        if (0 < angle < 90.0):
71            if vh > 0:
72                bearing = 90.0 - abs(angle)
73            if vh < 0:
74                bearing = 270.0 - abs(angle)
75        if (-90 < angle < 0):
76            if vh < 0:
77                bearing = 90.0 - abs(angle)
78            if vh > 0:
79                bearing = 270.0 - abs(angle)
80        if angle == 0:
81            bearing = 0.0
82               
83        model_time.append(t)       
84        stages.append(w)
85        elevations.append(z)  #Should be constant over time
86        momenta.append(m)
87        velocity.append(vel)
88        xmom.append(uh)
89        ymom.append(vh)
90        bearings.append(bearing)
91        depths.append(depth)
92
93        if w-z > max_depth:
94            max_depth = w-z
95        if m > max_momentum:           
96            max_momentum = m
97        if vel > max_velocity:
98            max_velocity = vel
99
100    #Plot only those gauges that have been inundated by more than a threshold
101    #if max_depth < 0.2:
102    #    print 'Skipping gauge %d' %k
103    #    continue
104
105    ion()
106    hold(False)
107
108    if elevations[0] <= 0:
109        plot(model_time, stages, '-b')
110    else:   
111        plot(model_time, stages, '-b',
112             model_time, elevations, '-k')
113    #axis([0, 100, z-0.01, z+0.005])
114       
115    #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
116    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
117    title(name)
118
119    title('%s (stage)' %name)
120    xlabel('time [s]')
121    ylabel('elevation [m]')   
122    legend(('Stage', 'Bed = %.1f' %elevations[0]),
123           shadow=True,
124           loc='upper right')
125
126    if store is True: savefig('Gauge_%s_stage' %myloc)
127
128    raw_input('Next')
129
130    """
131    #Momentum plot
132    ion()
133    hold(False)
134    plot(model_time, momenta, '-r')
135    title(name)
136
137    title('%s (momentum)' %name)
138    xlabel('time [s]')
139    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
140    #savefig('Gauge_%d_momentum' %k)
141    if store is True: savefig('Gauge_%s_momentum' %myloc)
142   
143    raw_input('Next')
144    """
145
146    """
147    #Bearing plot
148    ion()
149    hold(False)
150    ax = plot(model_time, bearings, '-b', model_time, due_west, '-.b',
151         model_time, due_east, '-.b')
152    title(name)
153    ax = axis([0, maxT, 0, 360])
154    text(maxT+tstep, 90, 'East')
155    text(maxT+tstep, 270, 'West')
156    #majorLocator = MultipleLocator(3600)
157    #print 'major', majorLocator[1]
158    #ax.xaxis.set_major_locator(majorLocator) #'yticklabels', range(30,390,30))
159    # set(labels,color='g',rotation=45)
160
161    title('%s (bearing)' %name)
162    xlabel('time [s]')
163    ylabel(' atan(vh/uh) [degrees from North]')   
164    #savefig('Gauge_%d_bearing' %k)
165    if store is True: savefig('Gauge_%s_bearing' %myloc)
166   
167    raw_input('Next')
168    """
169
170    """
171    #Speed plot
172    ion()
173    hold(False)
174    plot(model_time, velocity, '-r')
175    title(name)
176
177    title('%s (velocity)' %name)
178    xlabel('time [s]')
179    ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')   
180    #savefig('Gauge_%d_speed' %k)
181    if store is True: savefig('Gauge_%s_speed' %myloc)
182   
183    raw_input('Next')
184    """
185   
186
187    whichone = '_%s' %myloc
188    thisfile = 'island_gauge'+whichone+'.csv'
189    fid = open(thisfile, 'w')
190    s = 'Time, Depth, Momentum, Velocity \n'
191    fid.write(s)
192    for i_t, i_d, i_m, i_vel in zip(model_time, depths, momenta, velocity):
193        s = '%.2f, %.2f, %.2f, %.2f\n' %(i_t, i_d, i_m, i_vel)
194        fid.write(s)
195       
196show()
197
198
Note: See TracBrowser for help on using the repository browser.