[2404] | 1 | """Read in sww file, interpolate at specified locations and plot time series |
---|
| 2 | |
---|
| 3 | """ |
---|
| 4 | |
---|
| 5 | import project |
---|
| 6 | from pyvolution.util import file_function |
---|
| 7 | #from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn |
---|
| 8 | from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn |
---|
| 9 | from pylab import * |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | |
---|
| 13 | #swwfile = project.newoutputname + '.sww' |
---|
[2604] | 14 | swwfile = project.outputname |
---|
[2404] | 15 | |
---|
[2604] | 16 | # Causes a memory error |
---|
| 17 | #swwfile = 'O:/2/cit/inundation/Gippsland Lakes/120106/lakes_100_628759.sww' |
---|
| 18 | |
---|
| 19 | swwfile = 'O:/2/cit/inundation/Gippsland Lakes/lakes_100_314920.sww' |
---|
| 20 | |
---|
| 21 | |
---|
[2404] | 22 | #Time interval to plot |
---|
| 23 | tmin = 13000 |
---|
| 24 | tmax = 21000 |
---|
| 25 | |
---|
| 26 | def 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 | |
---|
| 47 | gauges, lines, locations = get_gauges_from_file(project.gauge_filename) |
---|
| 48 | |
---|
| 49 | #Read model output |
---|
[2604] | 50 | quantities = ['stage', 'elevation'] |
---|
[2404] | 51 | f = file_function(swwfile, |
---|
| 52 | quantities = quantities, |
---|
| 53 | interpolation_points = gauges, |
---|
| 54 | verbose = True, |
---|
[2604] | 55 | use_cache = False) |
---|
[2404] | 56 | |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | from math import sqrt |
---|
| 60 | N = len(gauges) |
---|
| 61 | for 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 |
---|
[2885] | 70 | for i, t in enumerate(f.get_time()): # T is a list of times |
---|
[2404] | 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 | |
---|
[2604] | 93 | plot(model_time, stages, '-b') |
---|
[2404] | 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 | |
---|
| 110 | show() |
---|
| 111 | |
---|