source: anuga_work/production/onslow_2006/get_timeseries.py @ 4504

Last change on this file since 4504 was 3535, checked in by duncan, 18 years ago

change imports so reflect the new structure

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