source: anuga_work/production/karratha_2005/process_buildings.py @ 5563

Last change on this file since 5563 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: 3.8 KB
Line 
1"""Read in sww file, interpolate at specified locations and plot time series
2Additionally, store augmented building database with max inundation depths
3"""
4
5import project
6from anuga.pyvolution.util import file_function
7from anuga.pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
8from pylab import *
9
10
11plot = False
12   
13swwfile = project.outputname + '.sww'
14
15#MOST Boundary directly north of Legendre island
16lat = project.north
17lon = (project.p3[1] + project.p4[1])/2 
18z, easting, northing  = redfearn(lat, lon)
19
20#Time interval to plot
21tmin = 13000
22tmax = 21000
23
24def get_gauges_from_file(filename):
25    fid = open(filename)
26    lines = fid.readlines()
27    fid.close()
28
29    gauges = []
30    for line in lines[1:]:
31        fields = line.split(',')
32        lon = float(fields[4])
33        lat = float(fields[5])
34
35        z, easting, northing  = redfearn(lat, lon)       
36        gauges.append([easting, northing])
37
38    #Return gauges and raw data for subsequent storage   
39    return gauges, lines
40
41gauges, buildings = get_gauges_from_file(project.gauge_filename)
42
43#Read model output
44quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
45f = file_function(swwfile,
46                  quantities = quantities,
47                  interpolation_points = gauges,
48                  verbose = True,
49                  use_cache = True)
50
51#New header in spreadsheet
52buildings[0] = buildings[0].strip() +\
53               ',MAX INUNDATION DEPTH (m)' +\
54               ',MAX MOMENTUM (m^2/s)\n'           
55               
56from math import sqrt
57N = len(gauges)
58for k, g in enumerate(gauges):
59    if k%((N+10)/10)==0:
60        print 'Doing row %d of %d' %(k, N)
61
62    model_time = []
63    stages = []
64    elevations = []
65    momenta = []
66
67    max_depth = 0
68    max_momentum = 0   
69    for i, t in enumerate(f.get_time()):
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    #Augment building data
88    buildings[k+1] = buildings[k+1].strip() +\
89                     ',%f' %max_depth +\
90                     ',%f\n' %max_momentum
91   
92    #print buildings[k]
93   
94    if plot is True:
95        #Plot only those gauges that have been inundated by more than a threshold
96        if max_depth < 0.2:
97            print 'Skipping gauge %d' %k
98            continue
99
100        #FIXME Reduce to tmin and tmax
101        #i0 = model_timefind
102   
103        ion()
104        hold(False)
105
106        if elevations[0] < -10:
107            plot(model_time, stages, '-b')
108        else:   
109            plot(model_time, stages, '-b',
110                 model_time, elevations, '-k')
111        name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
112        title(name)
113
114        title('%s (stage)' %name)
115        xlabel('time [s]')
116        ylabel('elevation [m]')   
117        legend(('Stage', 'Bed = %.1f' %elevations[0]),
118               shadow=True,
119               loc='upper right')
120        savefig('Gauge_%d_stage' %k)
121
122        raw_input('Next')
123
124
125        #Momentum plot
126        ion()
127        hold(False)
128        plot(model_time, momenta, '-r')
129        title(name)
130
131        title('%s (momentum)' %name)
132        xlabel('time [s]')
133        ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
134        savefig('Gauge_%d_momentum' %k)
135
136        raw_input('Next')
137   
138
139#Store new building file with max depths added
140
141FN = 'inundation_augmented_' + project.gauge_filename
142fid = open(FN, 'w')
143for line in buildings:
144    fid.write(line)
145fid.close()   
146
147
148if plot is True:   
149    show()
150
Note: See TracBrowser for help on using the repository browser.