source: branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py @ 7794

Last change on this file since 7794 was 7794, checked in by sexton, 12 years ago

updates to bunbury scripts

File size: 8.1 KB
Line 
1# ---------------------------------------------------------------------------
2# This script reads in both the master STS (output from process_gems_grids.py)
3# and the gems boundary order file (defined as gems_order in project.py)
4# to generate a boundary STS file which contains data just at the boundary
5# points for all timesteps, which is input into run_model.py in ANUGA.
6# Written by Nariman Habili and Leharne Fountain, 2010
7# ---------------------------------------------------------------------------
8
9import os
10import numpy
11from Scientific.IO.NetCDF import NetCDFFile
12from anuga.utilities.polygon import read_polygon
13import sys
14# sys.path.append('W:\\sandpits\\lfountain\\anuga_work\\production\\mandurah_storm_surge_2009\\gems_comparison')
15#sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
16import project
17from Scientific.IO.NetCDF import NetCDFFile
18from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference
19
20#-----------------
21# Directory setup
22#-------------------
23state = 'western_australia'
24scenario_folder = 'bunbury_storm_surge_scenario_2009'
25event = '20100527_gcom_12min' # Baseline TC Alby event
26
27ENV_INUNDATIONHOME = 'INUNDATIONHOME'
28home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
29# GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS')
30anuga_folder = os.path.join(home, state, scenario_folder, 'anuga')
31boundaries_folder = os.path.join(anuga_folder, 'boundaries')
32event_folder = os.path.join(boundaries_folder, event)
33
34#--------------------------------------------------------
35# Define boundary points to extract from master STS file
36#--------------------------------------------------------
37gems_boundary = numpy.array(read_polygon(project.gems_order))
38#gems_boundary = numpy.array(read_polygon(project.gauges))
39number_of_points = gems_boundary.shape[0]
40print 'hello', number_of_points
41print 'points', gems_boundary
42
43#---------------------------------------------------------
44# Read master sts file generated from GEMS ascii grids
45#---------------------------------------------------------
46print "Reading master NetCDF file"
47infile = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'rl') #rl for large netcdf files
48
49grid_size = infile.grid_size
50ncols = infile.ncols
51nrows = infile.nrows
52no_data = infile.no_data
53refzone = infile.zone
54x_origin = infile.xllcorner
55y_origin = infile.yllcorner
56x = infile.variables['x'][:]
57y = infile.variables['y'][:]
58time = infile.variables['time'][:]
59elevation = infile.variables['elevation'][:]
60stage_array = infile.variables['stage'][:]
61xmomentum_array = infile.variables['xmomentum'][:]
62ymomentum_array = infile.variables['ymomentum'][:]
63number_of_timesteps = stage_array.shape[0]
64
65infile.close()
66
67#--------------------------------------------
68# Creating index from boundary points
69#--------------------------------------------
70print "Creating index from boundary points" # finding location of boundary points in master file
71
72origin = numpy.array([x_origin[0], y_origin[0]])
73index_array_2D = ((gems_boundary - origin)/grid_size).astype(int)
74assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols # These test that the boundary points lie within the data area
75assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows
76
77index_array_1D = numpy.empty(index_array_2D.shape[0], int)
78index_array_1D = index_array_2D[:,1]*ncols + index_array_2D[:,0]
79
80#--------------------------------------------------
81# Extract values from master sts at index points
82#--------------------------------------------------
83print "Extract values from master sts at index points"
84
85x = x.take(index_array_1D)
86y = y.take(index_array_1D)
87elevation = elevation.take(index_array_1D)
88stage = numpy.empty([number_of_timesteps, number_of_points], float)
89xmomentum = numpy.empty([number_of_timesteps, number_of_points], float)
90ymomentum = numpy.empty([number_of_timesteps, number_of_points], float)
91
92for i, k in enumerate(index_array_1D):
93    stage[:,i] = stage_array[:,k] 
94    xmomentum[:,i] = xmomentum_array[:,k] 
95    ymomentum[:,i] = ymomentum_array[:,k] 
96
97#-----------------------------
98# Create the boundary STS file
99#-----------------------------
100print "Creating the boundary STS NetCDF file"
101
102fid = NetCDFFile(os.path.join(event_folder, event + '_boundary.sts'), 'w')
103fid.institution = 'Geoscience Australia'
104fid.description = "description"
105fid.starttime = 0.0
106fid.ncols = ncols
107fid.nrows = nrows
108fid.grid_size = grid_size
109fid.no_data = no_data
110fid.createDimension('number_of_points', number_of_points)
111fid.createDimension('number_of_timesteps', number_of_timesteps)
112fid.createDimension('numbers_in_range', 2)
113
114fid.createVariable('x', 'd', ('number_of_points',))
115fid.createVariable('y', 'd', ('number_of_points',))
116fid.createVariable('elevation', 'd', ('number_of_points',))
117fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
118fid.createVariable('time', 'i', ('number_of_timesteps',))
119fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
120fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
121fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
122fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
123fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
124fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
125
126fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
127fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
128fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
129fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
130fid.variables['elevation'][:] = elevation
131fid.variables['time'][:] = time
132
133s = fid.variables['stage']
134xm = fid.variables['xmomentum']
135ym = fid.variables['ymomentum']
136
137for i in xrange(number_of_timesteps):
138    s[i] = stage[i,:]
139    xm[i] = xmomentum[i,:]
140    ym[i] = ymomentum[i,:] 
141
142#---------------------------------------------------------------------------
143# Get timeseries values for wave height and components of momentum
144#---------------------------------------------------------------------------
145
146basename = 'sts_gauge_tide'
147for i in xrange(number_of_points):
148    out_file = os.path.join(event_folder,
149                            basename+'_'+str(i)+'.csv')
150    fid_sts = open(out_file, 'w')
151    fid_sts.write('time, stage, xmomentum, ymomentum, elevation \n')
152    for j in xrange(number_of_timesteps):
153        fid_sts.write('%.6f, %.6f, %.6f, %.6f, %.6f\n'
154                    % (time[j], stage[j,i], xmomentum[j,i], ymomentum[j,i], elevation[i]))
155
156    fid_sts.close() 
157
158origin = Geo_reference(refzone, x_origin, y_origin)
159geo_ref = write_NetCDF_georeference(origin, fid)
160
161fid.variables['x'][:] = x
162fid.variables['y'][:] = y 
163
164fid.close()
165
166sys.exit()
167
168print 'get timeseries'
169#-------------------------------------------------------------------------------
170# Get gauges (timeseries of index points)
171#-------------------------------------------------------------------------------
172
173basename = 'sts_gauge'
174quantity_names = ['stage', 'xmomentum', 'ymomentum']
175quantities = {}
176for i, name in enumerate(quantity_names):
177    quantities[name] = fid.variables[name][:]
178
179
180for j in range(len(x)):
181    index = j # permutation[j]
182    stage = quantities['stage'][:,j]
183    #print 'hello', j, stage
184    xmomentum = quantities['xmomentum'][:,j]
185    ymomentum = quantities['ymomentum'][:,j]
186
187    out_file = os.path.join(event_folder,
188                            basename+'_'+str(index)+'.csv')
189    fid_sts = open(out_file, 'w')
190    fid_sts.write('time, stage, xmomentum, ymomentum \n')
191
192    #-----------------------------------------------------------------------
193    # End of the get gauges
194    #-----------------------------------------------------------------------
195    for k in range(len(time)-1):
196        fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
197                      % (time[k], stage[k], xmomentum[k], ymomentum[k]))
198
199   
Note: See TracBrowser for help on using the repository browser.