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

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

updates to bunbury storm surge model

File size: 6.0 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')
15sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
16import project_GEMS_grid
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 = 'case_1' # 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_grid.gems_order))
38number_of_points = gems_boundary.shape[0]
39
40#---------------------------------------------------------
41# Read master sts file generated from GEMS ascii grids
42#---------------------------------------------------------
43print "Reading master NetCDF file"
44infile = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'r')
45
46grid_size = infile.grid_size
47ncols = infile.ncols
48nrows = infile.nrows
49no_data = infile.no_data
50refzone = infile.zone
51x_origin = infile.xllcorner
52y_origin = infile.yllcorner
53x = infile.variables['x'][:]
54y = infile.variables['y'][:]
55time = infile.variables['time'][:]
56elevation = infile.variables['elevation'][:]
57stage_array = infile.variables['stage'][:]
58xmomentum_array = infile.variables['xmomentum'][:]
59ymomentum_array = infile.variables['ymomentum'][:]
60number_of_timesteps = stage_array.shape[0]
61
62infile.close()
63
64#--------------------------------------------
65# Creating index from boundary points
66#--------------------------------------------
67print "Creating index from boundary points" # finding location of boundary points in master file
68
69origin = numpy.array([x_origin[0], y_origin[0]])
70index_array_2D = ((gems_boundary - origin)/grid_size).astype(int)
71assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols
72assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows
73
74index_array_1D = numpy.empty(index_array_2D.shape[0], int)
75index_array_1D = index_array_2D[:,1]*ncols + index_array_2D[:,0]
76
77#--------------------------------------------------
78# Extract values from master sts at index points
79#--------------------------------------------------
80print "Extract values from master sts at index points"
81
82x = x.take(index_array_1D)
83y = y.take(index_array_1D)
84elevation = elevation.take(index_array_1D)
85stage = numpy.empty([number_of_timesteps, number_of_points], float)
86xmomentum = numpy.empty([number_of_timesteps, number_of_points], float)
87ymomentum = numpy.empty([number_of_timesteps, number_of_points], float)
88
89for i, k in enumerate(index_array_1D):
90    stage[:,i] = stage_array[:,k] 
91    xmomentum[:,i] = xmomentum_array[:,k] 
92    ymomentum[:,i] = ymomentum_array[:,k] 
93
94#-----------------------------
95# Create the boundary STS file
96#-----------------------------
97print "Creating the boundary STS NetCDF file"
98
99fid = NetCDFFile(os.path.join(event_folder, event + '_boundary.sts'), 'w')
100fid.institution = 'Geoscience Australia'
101fid.description = "description"
102fid.starttime = 0.0
103fid.ncols = ncols
104fid.nrows = nrows
105fid.grid_size = grid_size
106fid.no_data = no_data
107fid.createDimension('number_of_points', number_of_points)
108fid.createDimension('number_of_timesteps', number_of_timesteps)
109fid.createDimension('numbers_in_range', 2)
110
111fid.createVariable('x', 'd', ('number_of_points',))
112fid.createVariable('y', 'd', ('number_of_points',))
113fid.createVariable('elevation', 'd', ('number_of_points',))
114fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
115fid.createVariable('time', 'i', ('number_of_timesteps',))
116fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
117fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
118fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
119fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
120fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
121fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))
122
123fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
124fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
125fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
126fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
127fid.variables['elevation'][:] = elevation
128fid.variables['time'][:] = time
129
130s = fid.variables['stage']
131xm = fid.variables['xmomentum']
132ym = fid.variables['ymomentum']
133
134for i in xrange(number_of_timesteps):
135    s[i] = stage[i,:]
136    xm[i] = xmomentum[i,:]
137    ym[i] = ymomentum[i,:]
138 
139origin = Geo_reference(refzone, x_origin, y_origin)
140geo_ref = write_NetCDF_georeference(origin, fid)
141
142fid.variables['x'][:] = x
143fid.variables['y'][:] = y
144
145fid.close()
Note: See TracBrowser for help on using the repository browser.