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

Last change on this file since 7688 was 7688, checked in by fountain, 14 years ago

updates to bunbury storm surge model

File size: 10.4 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_txt19780404.0*.gz'))
45grid_file_names.sort()
46
47# grid_file_names_temp = []
48# for i in range(0, 120, 4):
49    # grid_file_names_temp.append(grid_file_names[i])
50
51# grid_file_names = grid_file_names_temp
52   
53# number_of_timesteps = len(grid_file_names)
54
55print "Number of timesteps: ", number_of_timesteps
56
57start_time = 0      # all times in seconds
58timestep = 720     # all times in seconds
59end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
60refzone = 50        # UTM zone of model
61event_sts = os.path.join(event_folder, event)
62
63elevation = []
64stage = []
65speed = []
66theta = []
67
68elevation_file = open(elevation_file_name, 'rb')
69lines = elevation_file.readlines()
70elevation_file.close()
71
72# Strip off the ASCII header and also read ncols, nrows,
73# x_origin, y_origin, grid_size, no_data value, and read in elevation grid:
74for i, L in enumerate(lines):
75    if i == 0:
76        ncols = int(L.strip().split()[1])
77    if i == 1:
78        nrows = int(L.strip().split()[1])
79    if i == 2:
80        x_origin = int(L.strip().split()[1])
81    if i == 3:
82        y_origin = int(L.strip().split()[1])
83    if i == 4:
84        grid_size = int(L.strip().split()[1])
85    if i == 5:
86        no_data = int(float(L.strip().split()[1]))
87    if i > 5:
88        elevation+= L.strip().split()
89
90print 'Number or columns: ', ncols
91print 'Number of rows: ', nrows
92print 'X origin: ', x_origin
93print 'Y origin: ', y_origin
94print 'Grid size: ', grid_size
95print 'No data value: ', no_data
96       
97for f in grid_file_names:
98    print 'file: ', f
99   
100    gz_file = gzip.open(f, 'rb')
101    lines = gz_file.readlines()
102    gz_file.close()
103
104    for i, L in enumerate(lines):
105       
106        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
107            stage += [L.strip().split()]
108           
109        if i > (9 + nrows) and i < (10 + 2*nrows):
110            speed += [L.strip().split()]
111           
112        if i > (11 + 2*nrows):
113            theta += [L.strip().split()]
114
115#------------------------------------------------------
116# Create arrays of elevation, stage, speed and theta
117#-------------------------------------------------------
118print 'creating numpy arrays: elevation, stage, speed, theta'
119           
120elevation = numpy.array(elevation).astype('d')
121
122stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
123speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
124theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
125
126for i in range(number_of_timesteps):
127    ii = nrows*i
128    stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d'))
129    speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d'))
130    theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d'))
131
132stage = stage_array.ravel()
133speed = speed_array.ravel()
134theta = theta_array.ravel()
135   
136print 'Size of elevation array: ', elevation.size
137print 'Size of stage array: ', stage.size
138print 'Size of speed array: ', speed.size
139print 'Size of theta array: ', theta.size
140
141assert stage.size == speed.size == theta.size == ncols * nrows * number_of_timesteps
142assert stage.size == number_of_timesteps * elevation.size
143
144d = numpy.empty(elevation.size, dtype='d')
145depth = numpy.empty(stage.size, dtype='d')
146number_of_points = ncols * nrows
147
148# ---------------------------------------------
149# Create mask of no_data values across stage, speed and theta to ensure all three quantities
150# have no_data values in corresponding cells - this is a work-around until GEMS can produce a
151# a consistant dataset
152# ---------------------------------------------
153
154for i in xrange(number_of_timesteps):
155    j = number_of_points * i
156    k = j + number_of_points
157    # no_value_mask = (stage[j:k] < -9000) + (speed[j:k] < -9000) + (theta[j:k] < -9000)
158    no_value_index = numpy.where(((stage[j:k] < -100) + (speed[j:k] < 0) + (theta[j:k] < 0)) == True)[0]
159    depth[j:k] = stage[j:k] - elevation
160    numpy.put(stage[j:k], no_value_index, -9999)
161    numpy.put(speed[j:k], no_value_index, -9999)
162    numpy.put(theta[j:k], no_value_index, -9999)
163    numpy.put(depth[j:k], no_value_index, 0)
164
165# for i in xrange(number_of_timesteps):
166    # j = number_of_points * i
167    # k = j + number_of_points
168    # d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
169    # depth[j:k] = [0 if x<-9000 else x for x in d]   # sets stage to zero for no-data values of stage
170   
171    # 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
172    # for q in xrange(number_of_points):
173        # if elevation[q] >= 0 and stage[j+q]==0:
174            # depth[j+q] = 0
175   
176momentum = depth * speed
177
178xmomentum = numpy.empty(momentum.size, dtype='d')
179ymomentum = numpy.empty(momentum.size, dtype='d')
180
181print 'Calculating xmomentum and ymomentum'
182
183for i, t in enumerate(theta):
184   
185    mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
186    my = momentum[i]*math.cos(math.radians(t))
187
188    if t > 0.0 and t < 90.0:
189        assert mx > 0 and my > 0
190       
191    elif t > 90.0 and t < 180.0:
192        assert mx > 0 and my < 0
193       
194    elif t > 180.0 and t < 270.0:
195        assert mx < 0 and my < 0
196       
197    elif t > 270.0 and t < 360.0:
198        assert mx < 0 and my > 0
199       
200    elif t == 0.0 or t == 360.0:
201        assert my > 0
202       
203    elif t == 90.0:
204        assert mx > 0
205       
206    elif t == 180.0:
207        assert my < 0
208       
209    elif t == 270.0:
210        assert mx < 0
211       
212    elif t == -9999:
213        mx = 0
214        my = 0
215    else:
216        print "Unexpected value of theta"
217        exit()
218       
219    xmomentum[i] = mx
220    ymomentum[i] = my
221
222    assert math.sqrt(mx**2 + my**2) - momentum[i] < 1.e-06
223   
224x_origin_int = int(10000*x_origin)
225y_origin_int = int(10000*y_origin)
226grid_size_int = int(10000*grid_size)
227
228x = numpy.tile(numpy.arange(x_origin_int, (x_origin_int + ncols * grid_size_int), grid_size_int)/10000.0, nrows)
229y = numpy.repeat(numpy.arange(y_origin_int, (y_origin_int + nrows * grid_size_int), grid_size_int)/10000.0, ncols)
230
231assert x.size == y.size == number_of_points
232
233time = numpy.arange(start_time, end_time, timestep, dtype='i')
234
235assert time.size == number_of_timesteps
236assert momentum.size == stage.size
237
238# -----------------------------
239# Create the STS file
240# -----------------------------
241print "Creating the STS NetCDF file"
242
243fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w')
244fid.institution = 'Geoscience Australia'
245fid.description = "description"
246fid.starttime = 0.0
247fid.ncols = ncols
248fid.nrows = nrows
249fid.grid_size = grid_size
250fid.no_data = no_data
251fid.createDimension('number_of_points', number_of_points)
252fid.createDimension('number_of_timesteps', number_of_timesteps)
253fid.createDimension('numbers_in_range', 2)
254
255fid.createVariable('x', 'd', ('number_of_points',))
256fid.createVariable('y', 'd', ('number_of_points',))
257fid.createVariable('elevation', 'd', ('number_of_points',))
258fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
259fid.createVariable('time', 'i', ('number_of_timesteps',))
260fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
261fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
262fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
263fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
264fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
265fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
266
267fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
268fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
269fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
270fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
271fid.variables['elevation'][:] = elevation
272fid.variables['time'][:] = time
273
274s = fid.variables['stage']
275xm = fid.variables['xmomentum']
276ym = fid.variables['ymomentum']
277
278for i in xrange(number_of_timesteps):
279    ii = i*number_of_points
280    s[i] = stage[ii : ii + number_of_points]
281    xm[i] = xmomentum[ii : ii + number_of_points]
282    ym[i] = ymomentum[ii : ii + number_of_points]
283       
284origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings - it works for long and lat
285geo_ref = write_NetCDF_georeference(origin, fid)
286
287fid.variables['x'][:] = x - geo_ref.get_xllcorner()
288fid.variables['y'][:] = y - geo_ref.get_yllcorner()
289
290fid.close()
Note: See TracBrowser for help on using the repository browser.