source: branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py @ 7683

Last change on this file since 7683 was 7683, checked in by fountain, 14 years ago
File size: 9.7 KB
Line 
1# ---------------------------------------------------------------------------
2# This python script reads in GEMSURGE outputs and writes them into a master
3# STS file that contains all the data across the entire model domain at all
4# timesteps. It requires a vertically flipped elevation ASCII grid which is
5# merged with the GEMSURGE data to calculate the required quantities of
6# xmomentum and ymomentum.
7# Written by Nariman Habili and Leharne Fountain, 2010
8# ---------------------------------------------------------------------------
9
10import glob
11import os
12import gzip
13import pdb 
14import numpy
15import math
16from Scientific.IO.NetCDF import NetCDFFile
17from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference
18
19#-----------------
20# Directory setup
21#-------------------
22state = 'western_australia'
23scenario_folder = 'bunbury_storm_surge_scenario_2009'
24event = 'case_1'
25
26print "Processing event: ", event
27
28ENV_INUNDATIONHOME = 'INUNDATIONHOME'
29home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
30# GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS')
31anuga_folder = os.path.join(home, state, scenario_folder, 'anuga')
32boundaries_folder = os.path.join(anuga_folder, 'boundaries')
33event_folder = os.path.join(boundaries_folder, event)
34
35# input_dir = os.path.join(GEMS_folder, event)
36
37#-----------------------
38# Input files from GEMS
39#-----------------------
40elevation_file_name = os.path.join(event_folder, 'buntopog_20m_flip.asc')  # Name of the vertically flipped elevation grid
41
42print "Elevation file name: ", elevation_file_name
43
44grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780403.18*.gz'))
45grid_file_names.sort()
46number_of_timesteps = len(grid_file_names)
47
48print "Number of timesteps: ", number_of_timesteps
49
50start_time = 0      # all times in seconds
51timestep = 720     # all times in seconds
52end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
53refzone = 50        # UTM zone of model
54event_sts = os.path.join(event_folder, event)
55
56elevation = []
57stage = []
58speed = []
59theta = []
60
61elevation_file = open(elevation_file_name, 'rb')
62lines = elevation_file.readlines()
63elevation_file.close()
64
65# Strip off the ASCII header and also read ncols, nrows,
66# x_origin, y_origin, grid_size, no_data value, and read in elevation grid:
67for i, L in enumerate(lines):
68    if i == 0:
69        ncols = int(L.strip().split()[1])
70    if i == 1:
71        nrows = int(L.strip().split()[1])
72    if i == 2:
73        x_origin = int(L.strip().split()[1])
74    if i == 3:
75        y_origin = int(L.strip().split()[1])
76    if i == 4:
77        grid_size = int(L.strip().split()[1])
78    if i == 5:
79        no_data = int(float(L.strip().split()[1]))
80    if i > 5:
81        elevation+= L.strip().split()
82
83print 'Number or columns: ', ncols
84print 'Number of rows: ', nrows
85print 'X origin: ', x_origin
86print 'Y origin: ', y_origin
87print 'Grid size: ', grid_size
88print 'No data value: ', no_data
89       
90for f in grid_file_names:
91    print 'file: ', f
92   
93    gz_file = gzip.open(f, 'rb')
94    lines = gz_file.readlines()
95    gz_file.close()
96
97    for i, L in enumerate(lines):
98       
99        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
100            stage += [L.strip().split()]
101           
102        if i > (9 + nrows) and i < (10 + 2*nrows):
103            speed += [L.strip().split()]
104           
105        if i > (11 + 2*nrows):
106            theta += [L.strip().split()]
107
108#------------------------------------------------------
109# Create arrays of elevation, stage, speed and theta
110#-------------------------------------------------------
111print 'creating numpy arrays: elevation, stage, speed, theta'
112           
113elevation = numpy.array(elevation).astype('d')
114
115stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
116speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
117theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
118
119for i in range(number_of_timesteps):
120    ii = nrows*i
121    stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d'))
122    speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d'))
123    theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d'))
124
125stage = stage_array.ravel()
126speed = speed_array.ravel()
127theta = theta_array.ravel()
128   
129print 'Size of elevation array: ', elevation.size
130print 'Size of stage array: ', stage.size
131print 'Size of speed array: ', speed.size
132print 'Size of theta array: ', theta.size
133
134assert stage.size == speed.size == theta.size == ncols * nrows * number_of_timesteps
135assert stage.size == number_of_timesteps * elevation.size
136
137d = numpy.empty(elevation.size, dtype='d')
138depth = numpy.empty(stage.size, dtype='d')
139number_of_points = ncols * nrows
140
141for i in xrange(number_of_timesteps):
142    j = number_of_points * i
143    k = j + number_of_points
144    d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
145    depth[j:k] = [0 if x<-500 else x for x in d]   # sets stage to zero for no-data values of stage
146   
147    # depth[j:k] = [x if x>=0 and elevation>0 else 0 for x in d] # This is to correct for GEMS data which assumes stage=0 onshore when dry but should = elevation
148    # for q in xrange(number_of_points):
149        # if elevation[q] >= 0 and stage[j+q]==0:
150            # depth[j+q] = 0
151   
152momentum = depth * speed
153
154xmomentum = numpy.empty(momentum.size, dtype='d')
155ymomentum = numpy.empty(momentum.size, dtype='d')
156
157print 'Calculating xmomentum and ymomentum'
158
159for i, t in enumerate(theta):
160
161    if t > 0.0 and t < 90.0:
162        mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
163        my = momentum[i]*math.cos(math.radians(t))
164        assert mx >= 0 and my >= 0
165    if t > 90.0 and t < 180.0:
166        mx = momentum[i]*math.sin(math.radians(t))
167        my = momentum[i]*math.cos(math.radians(t))
168        assert mx >= 0 and my <= 0
169    if t > 180.0 and t < 270.0:
170        mx = momentum[i]*math.cos(math.radians(t))
171        my = momentum[i]*math.sin(math.radians(t))
172        assert mx <= 0 and my <= 0
173    if t > 270.0 and t < 360.0:
174        mx = momentum[i]*math.sin(math.radians(t))
175        my = momentum[i]*math.cos(math.radians(t))
176        assert mx <= 0 and my >= 0
177    if t == 0.0 or t == 360.0:
178        mx = 0
179        my = momentum[i]
180        assert my >= 0
181    if t == 90.0:
182        mx = momentum[i]
183        my = 0
184        assert mx >= 0
185    if t == 180.0:
186        mx = 0
187        my = -momentum[i]
188        assert my <= 0
189    if t == 270.0:
190        mx = -momentum[i]
191        my = 0
192        assert mx <= 0
193    if t == -9999:
194        mx = 0
195        my = 0
196       
197    xmomentum[i] = mx
198    ymomentum[i] = my
199
200    assert math.sqrt(mx**2 + my**2) - momentum[i] < 1.e-06
201   
202x_origin_int = int(10000*x_origin)
203y_origin_int = int(10000*y_origin)
204grid_size_int = int(10000*grid_size)
205
206x = numpy.tile(numpy.arange(x_origin_int, (x_origin_int + ncols * grid_size_int), grid_size_int)/10000.0, nrows)
207y = numpy.repeat(numpy.arange(y_origin_int, (y_origin_int + nrows * grid_size_int), grid_size_int)/10000.0, ncols)
208
209assert x.size == y.size == number_of_points
210
211time = numpy.arange(start_time, end_time, timestep, dtype='i')
212
213assert time.size == number_of_timesteps
214assert momentum.size == stage.size
215
216#-----------------------------
217# Create the STS file
218#-----------------------------
219print "Creating the STS NetCDF file"
220
221fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w')
222fid.institution = 'Geoscience Australia'
223fid.description = "description"
224fid.starttime = 0.0
225fid.ncols = ncols
226fid.nrows = nrows
227fid.grid_size = grid_size
228fid.no_data = no_data
229fid.createDimension('number_of_points', number_of_points)
230fid.createDimension('number_of_timesteps', number_of_timesteps)
231fid.createDimension('numbers_in_range', 2)
232
233fid.createVariable('x', 'd', ('number_of_points',))
234fid.createVariable('y', 'd', ('number_of_points',))
235fid.createVariable('elevation', 'd', ('number_of_points',))
236fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
237fid.createVariable('time', 'i', ('number_of_timesteps',))
238fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
239fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
240fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
241fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
242fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
243fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
244
245fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
246fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
247fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
248fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
249fid.variables['elevation'][:] = elevation
250fid.variables['time'][:] = time
251
252s = fid.variables['stage']
253xm = fid.variables['xmomentum']
254ym = fid.variables['ymomentum']
255
256for i in xrange(number_of_timesteps):
257    ii = i*number_of_points
258    s[i] = stage[ii : ii + number_of_points]
259    xm[i] = xmomentum[ii : ii + number_of_points]
260    ym[i] = ymomentum[ii : ii + number_of_points]
261       
262origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings - it works for long and lat
263geo_ref = write_NetCDF_georeference(origin, fid)
264
265fid.variables['x'][:] = x - geo_ref.get_xllcorner()
266fid.variables['y'][:] = y - geo_ref.get_yllcorner()
267
268fid.close()
Note: See TracBrowser for help on using the repository browser.