[1960] | 1 | """Read in sww file, read points along specified cross section and time and plot |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | import project |
---|
| 5 | from pyvolution.util import file_function |
---|
| 6 | from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | #swwfile = project.outputdir + project.basename + '.sww' |
---|
| 10 | #swwfile = project.outputname + '_0.0tide' + '.sww' |
---|
| 11 | swwfile = project.outputname + '.sww' |
---|
| 12 | swwfile = project.outputname + '_notsunami.sww' |
---|
| 13 | #swwfile = project.outputdir + 'karratha_100m_12m.sww' |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | time = 0 |
---|
| 17 | |
---|
| 18 | #E-W cross section |
---|
| 19 | start = [468500, 7715177] |
---|
| 20 | end = [469500, 7715177] |
---|
| 21 | step = 1 |
---|
| 22 | |
---|
| 23 | #Generate points (from west to east) |
---|
| 24 | points = [] |
---|
| 25 | x = start[0] |
---|
| 26 | y = start[1] |
---|
| 27 | while x <= end[0]: |
---|
| 28 | points.append( [x,y] ) |
---|
| 29 | x += step |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | #Read model output |
---|
| 33 | from caching import cache |
---|
| 34 | f = cache(file_function, swwfile, |
---|
| 35 | {'quantities': ['stage', 'elevation', 'xmomentum', 'ymomentum'], |
---|
| 36 | 'interpolation_points': points, |
---|
| 37 | 'verbose': True}, |
---|
| 38 | compression = False, |
---|
| 39 | dependencies = [swwfile], |
---|
| 40 | verbose = True) |
---|
| 41 | |
---|
[2885] | 42 | #print f |
---|
[1960] | 43 | |
---|
[2885] | 44 | T = f.get_time() |
---|
[1960] | 45 | from math import sqrt |
---|
[2885] | 46 | for time in [T[0], T[0] + 10]: |
---|
[1960] | 47 | |
---|
| 48 | abscissa = [] |
---|
| 49 | stages = [] |
---|
| 50 | elevations = [] |
---|
| 51 | momenta = [] |
---|
| 52 | |
---|
| 53 | for k, p in enumerate(points): |
---|
| 54 | # |
---|
| 55 | #print time, k, p |
---|
| 56 | w = f(time, point_id = k)[0] |
---|
| 57 | z = f(time, point_id = k)[1] |
---|
| 58 | uh = f(time, point_id = k)[2] |
---|
| 59 | vh = f(time, point_id = k)[3] |
---|
| 60 | |
---|
| 61 | abscissa.append( points[k][0] ) |
---|
| 62 | stages.append(w) |
---|
| 63 | elevations.append(z) |
---|
| 64 | momenta.append( sqrt(uh*uh + vh*vh) ) #Absolute momentum |
---|
| 65 | |
---|
| 66 | from pylab import * |
---|
| 67 | |
---|
| 68 | ion() |
---|
| 69 | hold(False) |
---|
| 70 | |
---|
| 71 | if elevations[0] < -10: |
---|
| 72 | plot(abscissa, stages, '-b') |
---|
| 73 | else: |
---|
| 74 | plot(abscissa, stages, '-b', |
---|
| 75 | abscissa, elevations, '-k') |
---|
| 76 | name = 'Cross section (time = %.1f)' %(time) |
---|
| 77 | title(name) |
---|
| 78 | |
---|
| 79 | title('%s (stage)' %name) |
---|
| 80 | xlabel('time [s]') |
---|
| 81 | ylabel('elevation [m]') |
---|
| 82 | legend(('Stage', 'Bed = %.1f' %elevations[0]), |
---|
| 83 | shadow=True, |
---|
| 84 | loc='upper right') |
---|
| 85 | savefig('Cross_section_%d' %int(time)) |
---|
| 86 | |
---|
| 87 | raw_input('Next') |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | #Momentum plot |
---|
| 91 | ion() |
---|
| 92 | hold(False) |
---|
| 93 | plot(abscissa, momenta, '-r') |
---|
| 94 | title(name) |
---|
| 95 | |
---|
| 96 | title('%s (momentum)' %name) |
---|
| 97 | xlabel('time [s]') |
---|
| 98 | ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]') |
---|
| 99 | savefig('Cross_section_%d_momentum' %int(time)) |
---|
| 100 | |
---|
| 101 | raw_input('Next') |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | show() |
---|
| 106 | |
---|