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 | |
---|
9 | import os |
---|
10 | import numpy |
---|
11 | from Scientific.IO.NetCDF import NetCDFFile |
---|
12 | from anuga.utilities.polygon import read_polygon |
---|
13 | import 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') |
---|
16 | import project_GEMS_grid |
---|
17 | from Scientific.IO.NetCDF import NetCDFFile |
---|
18 | from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference |
---|
19 | |
---|
20 | #----------------- |
---|
21 | # Directory setup |
---|
22 | #------------------- |
---|
23 | state = 'western_australia' |
---|
24 | scenario_folder = 'bunbury_storm_surge_scenario_2009' |
---|
25 | event = 'case_1' # Baseline TC Alby event |
---|
26 | |
---|
27 | ENV_INUNDATIONHOME = 'INUNDATIONHOME' |
---|
28 | home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder |
---|
29 | # GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS') |
---|
30 | anuga_folder = os.path.join(home, state, scenario_folder, 'anuga') |
---|
31 | boundaries_folder = os.path.join(anuga_folder, 'boundaries') |
---|
32 | event_folder = os.path.join(boundaries_folder, event) |
---|
33 | |
---|
34 | #-------------------------------------------------------- |
---|
35 | # Define boundary points to extract from master STS file |
---|
36 | #-------------------------------------------------------- |
---|
37 | gems_boundary = numpy.array(read_polygon(project_GEMS_grid.gems_order)) |
---|
38 | number_of_points = gems_boundary.shape[0] |
---|
39 | |
---|
40 | #--------------------------------------------------------- |
---|
41 | # Read master sts file generated from GEMS ascii grids |
---|
42 | #--------------------------------------------------------- |
---|
43 | print "Reading master NetCDF file" |
---|
44 | infile = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'r') |
---|
45 | |
---|
46 | grid_size = infile.grid_size |
---|
47 | ncols = infile.ncols |
---|
48 | nrows = infile.nrows |
---|
49 | no_data = infile.no_data |
---|
50 | refzone = infile.zone |
---|
51 | x_origin = infile.xllcorner |
---|
52 | y_origin = infile.yllcorner |
---|
53 | x = infile.variables['x'][:] |
---|
54 | y = infile.variables['y'][:] |
---|
55 | time = infile.variables['time'][:] |
---|
56 | elevation = infile.variables['elevation'][:] |
---|
57 | stage_array = infile.variables['stage'][:] |
---|
58 | xmomentum_array = infile.variables['xmomentum'][:] |
---|
59 | ymomentum_array = infile.variables['ymomentum'][:] |
---|
60 | number_of_timesteps = stage_array.shape[0] |
---|
61 | |
---|
62 | infile.close() |
---|
63 | |
---|
64 | #-------------------------------------------- |
---|
65 | # Creating index from boundary points |
---|
66 | #-------------------------------------------- |
---|
67 | print "Creating index from boundary points" # finding location of boundary points in master file |
---|
68 | |
---|
69 | origin = numpy.array([x_origin[0], y_origin[0]]) |
---|
70 | index_array_2D = ((gems_boundary - origin)/grid_size).astype(int) |
---|
71 | assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols |
---|
72 | assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows |
---|
73 | |
---|
74 | index_array_1D = numpy.empty(index_array_2D.shape[0], int) |
---|
75 | index_array_1D = index_array_2D[:,1]*ncols + index_array_2D[:,0] |
---|
76 | |
---|
77 | #-------------------------------------------------- |
---|
78 | # Extract values from master sts at index points |
---|
79 | #-------------------------------------------------- |
---|
80 | print "Extract values from master sts at index points" |
---|
81 | |
---|
82 | x = x.take(index_array_1D) |
---|
83 | y = y.take(index_array_1D) |
---|
84 | elevation = elevation.take(index_array_1D) |
---|
85 | stage = numpy.empty([number_of_timesteps, number_of_points], float) |
---|
86 | xmomentum = numpy.empty([number_of_timesteps, number_of_points], float) |
---|
87 | ymomentum = numpy.empty([number_of_timesteps, number_of_points], float) |
---|
88 | |
---|
89 | for 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 | #----------------------------- |
---|
97 | print "Creating the boundary STS NetCDF file" |
---|
98 | |
---|
99 | fid = NetCDFFile(os.path.join(event_folder, event + '_boundary.sts'), 'w') |
---|
100 | fid.institution = 'Geoscience Australia' |
---|
101 | fid.description = "description" |
---|
102 | fid.starttime = 0.0 |
---|
103 | fid.ncols = ncols |
---|
104 | fid.nrows = nrows |
---|
105 | fid.grid_size = grid_size |
---|
106 | fid.no_data = no_data |
---|
107 | fid.createDimension('number_of_points', number_of_points) |
---|
108 | fid.createDimension('number_of_timesteps', number_of_timesteps) |
---|
109 | fid.createDimension('numbers_in_range', 2) |
---|
110 | |
---|
111 | fid.createVariable('x', 'd', ('number_of_points',)) |
---|
112 | fid.createVariable('y', 'd', ('number_of_points',)) |
---|
113 | fid.createVariable('elevation', 'd', ('number_of_points',)) |
---|
114 | fid.createVariable('elevation_range', 'd', ('numbers_in_range',)) |
---|
115 | fid.createVariable('time', 'i', ('number_of_timesteps',)) |
---|
116 | fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points')) |
---|
117 | fid.createVariable('stage_range', 'd', ('numbers_in_range', )) |
---|
118 | fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points')) |
---|
119 | fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',)) |
---|
120 | fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points')) |
---|
121 | fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',)) |
---|
122 | |
---|
123 | fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
124 | fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
125 | fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
126 | fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
127 | fid.variables['elevation'][:] = elevation |
---|
128 | fid.variables['time'][:] = time |
---|
129 | |
---|
130 | s = fid.variables['stage'] |
---|
131 | xm = fid.variables['xmomentum'] |
---|
132 | ym = fid.variables['ymomentum'] |
---|
133 | |
---|
134 | for i in xrange(number_of_timesteps): |
---|
135 | s[i] = stage[i,:] |
---|
136 | xm[i] = xmomentum[i,:] |
---|
137 | ym[i] = ymomentum[i,:] |
---|
138 | |
---|
139 | origin = Geo_reference(refzone, x_origin, y_origin) |
---|
140 | geo_ref = write_NetCDF_georeference(origin, fid) |
---|
141 | |
---|
142 | fid.variables['x'][:] = x |
---|
143 | fid.variables['y'][:] = y |
---|
144 | |
---|
145 | fid.close() |
---|