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

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

updates to bunbury storm surge model

File size: 9.5 KB
RevLine 
[7678]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'
[7681]24event = 'case_1'
[7678]25
[7681]26print "Processing event: ", event
27
[7678]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#-----------------------
[7681]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*.gz'))
[7678]45grid_file_names.sort()
46number_of_timesteps = len(grid_file_names)
[7681]47grid_size = 20  # cellsize from the GEMSURGE output grids
48start_time = 0      # all times in seconds
49end_time = 61200    # end_time = start_time + timestep*number of grid_file_names. all times in seconds
50timestep = 720     # all times in seconds
51refzone = 50        # UTM zone of model
52x_origin = 368600 #364221.03   # xllcorner from GEMSURGE output grids
53y_origin = 6304600 #6371062.71  # yllcorner from GEMSURGE output grids
[7678]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,
[7681]66# no_data value, and read in elevation grid:
[7678]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])
[7681]72    if i == 5:
73        no_data = int(float(L.strip().split()[1]))
[7678]74    if i > 5:
75        elevation+= L.strip().split()
76
77print 'Number or columns: ', ncols
78print 'Number of rows: ', nrows
[7681]79       
[7678]80for f in grid_file_names:
81    print 'file: ', f
82   
83    gz_file = gzip.open(f, 'rb')
84    lines = gz_file.readlines()
85    gz_file.close()
86
87    for i, L in enumerate(lines):
88       
[7681]89        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
[7678]90            stage += [L.strip().split()]
91           
92        if i > (9 + nrows) and i < (10 + 2*nrows):
93            speed += [L.strip().split()]
94           
95        if i > (11 + 2*nrows):
96            theta += [L.strip().split()]
97
98#------------------------------------------------------
99# Create arrays of elevation, stage, speed and theta
100#-------------------------------------------------------
101print 'creating numpy arrays: elevation, stage, speed, theta'
102           
103elevation = numpy.array(elevation).astype('d')
104
105stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
106speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
107theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
108
109for i in range(number_of_timesteps):
110    ii = nrows*i
111    stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d'))
112    speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d'))
113    theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d'))
114
115stage = stage_array.ravel()
116speed = speed_array.ravel()
117theta = theta_array.ravel()
118   
119print 'Size of elevation array: ', elevation.size
120print 'Size of stage array: ', stage.size
121print 'Size of speed array: ', speed.size
122print 'Size of theta array: ', theta.size
123
124assert stage.size == speed.size == theta.size == ncols * nrows * number_of_timesteps
125assert stage.size == number_of_timesteps * elevation.size
126
[7681]127d = numpy.empty(elevation.size, dtype='d')
[7678]128depth = numpy.empty(stage.size, dtype='d')
129number_of_points = ncols * nrows
130
131for i in xrange(number_of_timesteps):
132    j = number_of_points * i
133    k = j + number_of_points
[7681]134    d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
135    depth[j:k] = [0 if x<-9900 else x for x in d]   # sets stage to zero for no-data values of stage
136   
137    # 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
138    # for q in xrange(number_of_points):
139        # if elevation[q] >= 0 and stage[j+q]==0:
140            # depth[j+q] = 0
141   
[7678]142momentum = depth * speed
143
144xmomentum = numpy.empty(momentum.size, dtype='d')
145ymomentum = numpy.empty(momentum.size, dtype='d')
146
147print 'Calculating xmomentum and ymomentum'
148
149for i, t in enumerate(theta):
150
151    if t > 0.0 and t < 90.0:
152        mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
153        my = momentum[i]*math.cos(math.radians(t))
154        assert mx > 0 and my > 0
155    if t > 90.0 and t < 180.0:
[7681]156        mx = momentum[i]*math.sin(math.radians(t))
157        my = momentum[i]*math.cos(math.radians(t))
[7678]158        assert mx > 0 and my < 0
159    if t > 180.0 and t < 270.0:
[7681]160        mx = momentum[i]*math.cos(math.radians(t))
161        my = momentum[i]*math.sin(math.radians(t))
[7678]162        assert mx < 0 and my < 0
163    if t > 270.0 and t < 360.0:
[7681]164        mx = momentum[i]*math.sin(math.radians(t))
165        my = momentum[i]*math.cos(math.radians(t))
[7678]166        assert mx < 0 and my > 0
167    if t == 0.0 or t == 360.0:
168        mx = 0
169        my = momentum[i]
[7681]170        assert my >= 0
[7678]171    if t == 90.0:
172        mx = momentum[i]
173        my = 0
[7681]174        assert mx >= 0
[7678]175    if t == 180.0:
176        mx = 0
[7681]177        my = -momentum[i]
178        assert my <= 0
[7678]179    if t == 270.0:
[7681]180        mx = -momentum[i]
[7678]181        my = 0
[7681]182        assert mx <= 0
[7678]183       
184    xmomentum[i] = mx
185    ymomentum[i] = my
186
[7681]187    assert math.sqrt(mx**2 + my**2) - momentum[i] < 1.e-06
[7678]188   
189x_origin_int = int(10000*x_origin)
190y_origin_int = int(10000*y_origin)
191grid_size_int = int(10000*grid_size)
192
193x = numpy.tile(numpy.arange(x_origin_int, (x_origin_int + ncols * grid_size_int), grid_size_int)/10000.0, nrows)
194y = numpy.repeat(numpy.arange(y_origin_int, (y_origin_int + nrows * grid_size_int), grid_size_int)/10000.0, ncols)
195
196assert x.size == y.size == number_of_points
197
198time = numpy.arange(start_time, end_time, timestep, dtype='i')
199
200assert time.size == number_of_timesteps
201assert momentum.size == stage.size
202
203#-----------------------------
204# Create the STS file
205#-----------------------------
206print "Creating the STS NetCDF file"
207
[7681]208fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w')
[7678]209fid.institution = 'Geoscience Australia'
210fid.description = "description"
211fid.starttime = 0.0
212fid.ncols = ncols
213fid.nrows = nrows
214fid.grid_size = grid_size
215fid.no_data = no_data
216fid.createDimension('number_of_points', number_of_points)
217fid.createDimension('number_of_timesteps', number_of_timesteps)
218fid.createDimension('numbers_in_range', 2)
219
220fid.createVariable('x', 'd', ('number_of_points',))
221fid.createVariable('y', 'd', ('number_of_points',))
222fid.createVariable('elevation', 'd', ('number_of_points',))
223fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
224fid.createVariable('time', 'i', ('number_of_timesteps',))
225fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
226fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
227fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
228fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
229fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
230fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
231
232fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
233fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
234fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
235fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
236fid.variables['elevation'][:] = elevation
237fid.variables['time'][:] = time
238
239s = fid.variables['stage']
240xm = fid.variables['xmomentum']
241ym = fid.variables['ymomentum']
242
243for i in xrange(number_of_timesteps):
244    ii = i*number_of_points
245    s[i] = stage[ii : ii + number_of_points]
246    xm[i] = xmomentum[ii : ii + number_of_points]
247    ym[i] = ymomentum[ii : ii + number_of_points]
248       
249origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings - it works for long and lat
250geo_ref = write_NetCDF_georeference(origin, fid)
251
252fid.variables['x'][:] = x - geo_ref.get_xllcorner()
253fid.variables['y'][:] = y - geo_ref.get_yllcorner()
254
255fid.close()
Note: See TracBrowser for help on using the repository browser.