Changeset 2069
- Timestamp:
- Nov 24, 2005, 6:12:22 PM (18 years ago)
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pyvolution/data_manager.py
r2064 r2069 362 362 from os import stat 363 363 364 minimum_allowed_depth = 0.001365 #minimum_allowed_depth = 0.0 #FIXME pass in or read from domain364 #minimum_allowed_depth = 0.001 365 minimum_allowed_depth = 0.0 #FIXME pass in or read from domain 366 366 from Numeric import choose 367 367 -
production/karratha_2005/get_timeseries.py
r1987 r2069 1 """Read in sww file, interpolate at specified locations and plot 1 """Read in sww file, interpolate at specified locations and plot time series 2 Additionally, store augmented building database with max inundation depths 2 3 """ 3 4 … … 6 7 from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn 7 8 from pylab import * 8 9 9 10 10 11 #swwfile = project.outputdir + project.basename + '.sww' 12 #swwfile = project.outputname + '_0.0tide' + '.sww' 11 plot = False 12 13 13 swwfile = project.outputname + '.sww' 14 #swwfile = project.outputname + '_notsunami.sww'15 #swwfile = project.outputdir + 'karratha_100m_12m.sww'16 14 17 15 #MOST Boundary directly north of Legendre island … … 32 30 for line in lines[1:]: 33 31 fields = line.split(',') 34 lon = float(fields[ 5])35 lat = float(fields[ 6])32 lon = float(fields[4]) 33 lat = float(fields[5]) 36 34 37 35 z, easting, northing = redfearn(lat, lon) 38 36 gauges.append([easting, northing]) 39 37 40 return gauges 38 #Return gauges and raw data for subsequent storage 39 return gauges, lines 41 40 42 gauges = get_gauges_from_file(project.gauge_filename) 43 44 45 #gauges = gauges[-10:] 46 #print gauges 47 #import sys; sys.exit() 48 49 #From Neil 50 #gauges = [[easting, northing - 5], 51 # [easting, northing - 10], 52 # [easting, northing - 50], 53 # [easting, northing - 100], 54 # [470882, 7717952], 55 # [469390, 7715426], 56 # [469214, 7714979], 57 # [468899, 7715177], 58 # [469038, 7715725], 59 # [470285, 7717470]] 60 41 gauges, buildings = get_gauges_from_file(project.gauge_filename) 61 42 62 43 #Read model output … … 68 49 use_cache = True) 69 50 70 print f71 #model_time = f.T72 51 73 52 buildings[0] = buildings[0].strip() +\ 53 ',MAX INUNDATION DEPTH (m)' +\ 54 ',MAX MOMENTUM (m^2/s)\n' 55 74 56 from math import sqrt 57 N = len(gauges) 75 58 for k, g in enumerate(gauges): 59 if k%((N+10)/10)==0: 60 print 'Doing row %d of %d' %(k, N) 76 61 77 62 model_time = [] … … 80 65 momenta = [] 81 66 67 max_depth = 0 68 max_momentum = 0 69 for i, t in enumerate(f.T): 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] 82 74 83 max_depth = 0 84 for i, t in enumerate(f.T): 85 if tmin < t < tmax: 86 w = f(t, point_id = k)[0] 87 z = f(t, point_id = k)[1] 88 uh = f(t, point_id = k)[2] 89 vh = f(t, point_id = k)[3] 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) 90 81 91 model_time.append(t)92 stages.append(w)93 elevations.append(z) #Should be constant over time94 m omenta.append( sqrt(uh*uh + vh*vh) ) #Absolute momentum82 if w-z > max_depth: 83 max_depth = w-z 84 if m > max_momentum: 85 max_momentum = m 95 86 96 if w-z > max_depth: 97 max_depth = w-z 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 98 99 99 100 #FIXME Reduce to tmin and tmax 101 #i0 = model_timefind 102 103 ion() 104 hold(False) 100 105 101 #Plot only those gauges that have been inundated by more than a threshold102 if max_depth < 0.2:103 print 'Skipping gauge %d' %k104 continue105 106 ion()107 hold(False)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) 108 113 109 if elevations[0] < -10:110 plot(model_time, stages, '-b')111 else:112 plot(model_time, stages, '-b',113 model_time, elevations, '-k')114 name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])115 title(name)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) 116 121 117 title('%s (stage)' %name) 118 xlabel('time [s]') 119 ylabel('elevation [m]') 120 legend(('Stage', 'Bed = %.1f' %elevations[0]), 121 shadow=True, 122 loc='upper right') 123 savefig('Gauge_%d_stage' %k) 124 125 raw_input('Next') 122 raw_input('Next') 126 123 127 124 128 #Momentum plot129 ion()130 hold(False)131 plot(model_time, momenta, '-r')132 title(name)125 #Momentum plot 126 ion() 127 hold(False) 128 plot(model_time, momenta, '-r') 129 title(name) 133 130 134 title('%s (momentum)' %name)135 xlabel('time [s]')136 ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')137 savefig('Gauge_%d_momentum' %k)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) 138 135 139 raw_input('Next')136 raw_input('Next') 140 137 141 138 142 143 show() 139 #Store new building file with mak depths added 144 140 141 FN = 'inundation_augmented_' + project.gauge_filename 142 fid = open(FN, 'w') 143 for line in buildings: 144 fid.write(line) 145 fid.close() 146 147 148 if plot is True: 149 show() 150 -
production/karratha_2005/project.py
r2017 r2069 31 31 outputname = outputdir + basename #Used by post processing 32 32 33 gauge_filename = 'all_bld_ind.csv' 33 #gauge_filename = 'all_bld_ind.csv' 34 gauge_filename = 'karratha_buildings_23112005.csv' 34 35 35 36
Note: See TracChangeset
for help on using the changeset viewer.