source: production/gippsland_2005/get_timeseries.py @ 2604

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

check in unknown changes

File size: 2.8 KB
Line 
1"""Read in sww file, interpolate at specified locations and plot time series
2
3"""
4
5import project
6from pyvolution.util import file_function
7#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
8from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
9from pylab import *
10
11
12   
13#swwfile = project.newoutputname + '.sww'
14swwfile = project.outputname
15
16# Causes a memory error
17#swwfile = 'O:/2/cit/inundation/Gippsland Lakes/120106/lakes_100_628759.sww'
18
19swwfile = 'O:/2/cit/inundation/Gippsland Lakes/lakes_100_314920.sww'
20
21
22#Time interval to plot
23tmin = 13000
24tmax = 21000
25
26def get_gauges_from_file(filename):
27    fid = open(filename)
28    lines = fid.readlines()
29    fid.close()
30   
31    gauges = []
32    gaugelocation = []
33    for line in lines[1:]:
34        fields = line.split(',')
35        # my gauge file set up as locationname, easting, northing
36        location = fields[0]
37        easting = float(fields[1])
38        northing = float(fields[2])
39        #z, easting, northing  = redfearn(lat, lon)
40        gauges.append([easting, northing])
41        gaugelocation.append(location)
42       
43    #Return gauges and raw data for subsequent storage
44    #return gauges, linesfs
45    return gauges, lines, gaugelocation
46
47gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
48
49#Read model output
50quantities = ['stage', 'elevation']
51f = file_function(swwfile,
52                  quantities = quantities,
53                  interpolation_points = gauges,
54                  verbose = True,
55                  use_cache = False)
56
57
58               
59from math import sqrt
60N = len(gauges)
61for k, g in enumerate(gauges):
62    if k%((N+10)/10)==0: # diagnostics - print 10 lines
63        print 'Doing row %d of %d' %(k, N)
64
65    model_time = []
66    stages = []
67    elevations = []
68
69    max_depth = 0
70    for i, t in enumerate(f.T): # T is a list of times
71        #if tmin < t < tmax:
72        w = f(t, point_id = k)[0]
73        z = f(t, point_id = k)[1]
74        #myloc = locations[k]
75       
76        model_time.append(t)       
77        stages.append(w)
78        elevations.append(z)  #Should be constant over time
79
80        if w-z > max_depth:
81            max_depth = w-z
82
83   
84
85    #Plot only those gauges that have been inundated by more than a threshold
86    #if max_depth < 0.2:
87    #    print 'Skipping gauge %d' %k
88    #    continue
89
90    ion()
91    hold(False)
92
93    plot(model_time, stages, '-b')     
94    name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
95    #name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
96    title(name)
97
98    title('%s (stage)' %name)
99    xlabel('time [s]')
100    ylabel('elevation [m]')   
101    legend(('Stage', 'Bed = %.1f' %elevations[0]),
102           shadow=True,
103           loc='upper right')
104    savefig('Gauge_%d_stage' %k) # savefig('Gauge_%s_stage' %myloc)
105
106    raw_input('Next')
107
108   
109
110show()
111
Note: See TracBrowser for help on using the repository browser.