[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 | |
---|
| 10 | import glob |
---|
| 11 | import os |
---|
| 12 | import gzip |
---|
| 13 | import pdb |
---|
| 14 | import numpy |
---|
| 15 | import math |
---|
| 16 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 17 | from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference |
---|
| 18 | |
---|
| 19 | #----------------- |
---|
| 20 | # Directory setup |
---|
| 21 | #------------------- |
---|
| 22 | state = 'western_australia' |
---|
| 23 | scenario_folder = 'bunbury_storm_surge_scenario_2009' |
---|
[7794] | 24 | event = '20100527_gcom_12min' |
---|
[7678] | 25 | |
---|
[7681] | 26 | print "Processing event: ", event |
---|
| 27 | |
---|
[7678] | 28 | ENV_INUNDATIONHOME = 'INUNDATIONHOME' |
---|
| 29 | home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder |
---|
| 30 | # GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS') |
---|
| 31 | anuga_folder = os.path.join(home, state, scenario_folder, 'anuga') |
---|
| 32 | boundaries_folder = os.path.join(anuga_folder, 'boundaries') |
---|
| 33 | event_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] | 40 | elevation_file_name = os.path.join(event_folder, 'buntopog_20m_flip.asc') # Name of the vertically flipped elevation grid |
---|
| 41 | |
---|
| 42 | print "Elevation file name: ", elevation_file_name |
---|
| 43 | |
---|
[7692] | 44 | grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz')) |
---|
[7678] | 45 | grid_file_names.sort() |
---|
[7683] | 46 | |
---|
[7794] | 47 | #grid_file_names_temp = [] |
---|
| 48 | # this was to subsample the gz files |
---|
| 49 | #for i in range(1, 32, 1): |
---|
| 50 | # grid_file_names_temp.append(grid_file_names[i]) |
---|
[7688] | 51 | |
---|
[7794] | 52 | #grid_file_names = grid_file_names_temp |
---|
[7688] | 53 | |
---|
[7692] | 54 | number_of_timesteps = len(grid_file_names) |
---|
[7688] | 55 | |
---|
[7683] | 56 | print "Number of timesteps: ", number_of_timesteps |
---|
| 57 | |
---|
[7681] | 58 | start_time = 0 # all times in seconds |
---|
[7794] | 59 | timestep = 720 # all times in seconds |
---|
[7683] | 60 | end_time = start_time + (timestep*number_of_timesteps) # all times in seconds |
---|
[7681] | 61 | refzone = 50 # UTM zone of model |
---|
[7678] | 62 | event_sts = os.path.join(event_folder, event) |
---|
| 63 | |
---|
| 64 | elevation = [] |
---|
| 65 | |
---|
| 66 | elevation_file = open(elevation_file_name, 'rb') |
---|
| 67 | lines = elevation_file.readlines() |
---|
| 68 | elevation_file.close() |
---|
| 69 | |
---|
| 70 | # Strip off the ASCII header and also read ncols, nrows, |
---|
[7683] | 71 | # x_origin, y_origin, grid_size, no_data value, and read in elevation grid: |
---|
[7678] | 72 | for i, L in enumerate(lines): |
---|
| 73 | if i == 0: |
---|
| 74 | ncols = int(L.strip().split()[1]) |
---|
| 75 | if i == 1: |
---|
| 76 | nrows = int(L.strip().split()[1]) |
---|
[7683] | 77 | if i == 2: |
---|
| 78 | x_origin = int(L.strip().split()[1]) |
---|
| 79 | if i == 3: |
---|
| 80 | y_origin = int(L.strip().split()[1]) |
---|
| 81 | if i == 4: |
---|
| 82 | grid_size = int(L.strip().split()[1]) |
---|
[7681] | 83 | if i == 5: |
---|
| 84 | no_data = int(float(L.strip().split()[1])) |
---|
[7678] | 85 | if i > 5: |
---|
| 86 | elevation+= L.strip().split() |
---|
| 87 | |
---|
| 88 | print 'Number or columns: ', ncols |
---|
| 89 | print 'Number of rows: ', nrows |
---|
[7683] | 90 | print 'X origin: ', x_origin |
---|
| 91 | print 'Y origin: ', y_origin |
---|
| 92 | print 'Grid size: ', grid_size |
---|
| 93 | print 'No data value: ', no_data |
---|
[7794] | 94 | |
---|
[7723] | 95 | #------------------------------------------------------ |
---|
[7794] | 96 | # Create arrays of elevation, depth, current_x and current_y |
---|
[7723] | 97 | #------------------------------------------------------- |
---|
[7794] | 98 | print 'creating numpy arrays: elevation, depth, current_x, current_y' |
---|
[7692] | 99 | |
---|
[7794] | 100 | depth = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float) |
---|
| 101 | current_x = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float) |
---|
| 102 | current_y = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float) |
---|
| 103 | |
---|
[7692] | 104 | for j, f in enumerate(grid_file_names): |
---|
[7678] | 105 | print 'file: ', f |
---|
| 106 | |
---|
| 107 | gz_file = gzip.open(f, 'rb') |
---|
| 108 | |
---|
[7794] | 109 | d = [] |
---|
| 110 | cx = [] |
---|
| 111 | cy = [] |
---|
| 112 | |
---|
[7692] | 113 | for i, L in enumerate(gz_file): |
---|
[7678] | 114 | |
---|
[7794] | 115 | #if i > 7 and i < (8 + nrows): # 7 refers to the number of rows of header info that we skip - always CHECK format |
---|
| 116 | # print 'first block' |
---|
| 117 | |
---|
| 118 | #if i > (9 + nrows) and i < (10 + 2*nrows): |
---|
| 119 | # print 'second block' |
---|
| 120 | |
---|
| 121 | #if i > (11 + 2*nrows): |
---|
| 122 | # print 'third block' |
---|
| 123 | |
---|
| 124 | if i > 13 + 3*nrows and i < (14 + 4*nrows): |
---|
| 125 | d += L.strip().split() |
---|
[7678] | 126 | |
---|
[7794] | 127 | if i > (15 + 4*nrows) and i < (16 + 5*nrows): |
---|
| 128 | cx += L.strip().split() |
---|
[7678] | 129 | |
---|
[7794] | 130 | if i > (17 + 5*nrows): |
---|
| 131 | cy += L.strip().split() |
---|
[7678] | 132 | |
---|
[7794] | 133 | depth[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(d).astype('d')) |
---|
| 134 | current_x[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(cx).astype('d')) |
---|
| 135 | current_y[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(cy).astype('d')) |
---|
| 136 | |
---|
| 137 | gz_file.close() |
---|
| 138 | |
---|
[7678] | 139 | elevation = numpy.array(elevation).astype('d') |
---|
| 140 | |
---|
| 141 | print 'Size of elevation array: ', elevation.size |
---|
[7794] | 142 | print 'Size of depth array: ', depth.size |
---|
| 143 | print 'Size of current x array: ', current_x.size |
---|
| 144 | print 'Size of current y array: ', current_y.size |
---|
[7678] | 145 | |
---|
[7794] | 146 | assert depth.size == current_x.size == current_y.size == ncols * nrows * number_of_timesteps |
---|
| 147 | assert depth.size == number_of_timesteps * elevation.size |
---|
[7678] | 148 | |
---|
[7794] | 149 | stage = numpy.empty(depth.size, dtype='d') |
---|
[7678] | 150 | number_of_points = ncols * nrows |
---|
| 151 | |
---|
[7688] | 152 | # --------------------------------------------- |
---|
[7794] | 153 | # Create mask of no_data values across depth, current_x and current_y to ensure all three quantities |
---|
[7688] | 154 | # have no_data values in corresponding cells - this is a work-around until GEMS can produce a |
---|
| 155 | # a consistant dataset |
---|
| 156 | # --------------------------------------------- |
---|
[7794] | 157 | |
---|
| 158 | no_value_index = numpy.where(((depth < -9000) + (current_x < -9000) + (current_y < -9000)) == True)[0] |
---|
| 159 | |
---|
[7692] | 160 | numpy.put(stage, no_value_index, -9999) |
---|
[7794] | 161 | numpy.put(current_x, no_value_index, -9999) |
---|
| 162 | numpy.put(current_y, no_value_index, -9999) |
---|
[7723] | 163 | numpy.put(depth, no_value_index, 0) |
---|
[7681] | 164 | |
---|
[7794] | 165 | # Taking absolute value is to account for -ve depths obtained when depth-elevation |
---|
[7692] | 166 | # is slightly -ve - why I don't know, possbly a rounding error? |
---|
[7794] | 167 | #momentum = numpy.absolute(depth * current_x) |
---|
[7678] | 168 | |
---|
[7794] | 169 | print 'Calculating stage, xmomentum and ymomentum' |
---|
[7678] | 170 | |
---|
[7794] | 171 | #stage = depth - numpy.tile(numpy.absolute(elevation), number_of_timesteps) |
---|
| 172 | stage = depth + numpy.tile(elevation, number_of_timesteps) |
---|
| 173 | xmomentum = current_x*depth #momentum*numpy.sin(numpy.radians(current_y)) |
---|
| 174 | ymomentum = current_y*depth #momentum*numpy.cos(numpy.radians(current_y)) |
---|
| 175 | |
---|
| 176 | numpy.put(xmomentum, no_value_index, 0) |
---|
| 177 | numpy.put(ymomentum, no_value_index, 0) |
---|
| 178 | |
---|
| 179 | #assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0 |
---|
[7678] | 180 | |
---|
| 181 | x_origin_int = int(10000*x_origin) |
---|
| 182 | y_origin_int = int(10000*y_origin) |
---|
| 183 | grid_size_int = int(10000*grid_size) |
---|
| 184 | |
---|
| 185 | x = numpy.tile(numpy.arange(x_origin_int, (x_origin_int + ncols * grid_size_int), grid_size_int)/10000.0, nrows) |
---|
| 186 | y = numpy.repeat(numpy.arange(y_origin_int, (y_origin_int + nrows * grid_size_int), grid_size_int)/10000.0, ncols) |
---|
| 187 | |
---|
| 188 | assert x.size == y.size == number_of_points |
---|
| 189 | |
---|
| 190 | time = numpy.arange(start_time, end_time, timestep, dtype='i') |
---|
| 191 | |
---|
| 192 | assert time.size == number_of_timesteps |
---|
[7794] | 193 | assert xmomentum.size == depth.size == ymomentum.size |
---|
[7678] | 194 | |
---|
[7688] | 195 | # ----------------------------- |
---|
[7678] | 196 | # Create the STS file |
---|
[7688] | 197 | # ----------------------------- |
---|
[7678] | 198 | print "Creating the STS NetCDF file" |
---|
| 199 | |
---|
[7794] | 200 | #for j in range(2): |
---|
| 201 | |
---|
[7723] | 202 | fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'wl') |
---|
[7678] | 203 | fid.institution = 'Geoscience Australia' |
---|
| 204 | fid.description = "description" |
---|
| 205 | fid.starttime = 0.0 |
---|
| 206 | fid.ncols = ncols |
---|
| 207 | fid.nrows = nrows |
---|
| 208 | fid.grid_size = grid_size |
---|
| 209 | fid.no_data = no_data |
---|
| 210 | fid.createDimension('number_of_points', number_of_points) |
---|
| 211 | fid.createDimension('number_of_timesteps', number_of_timesteps) |
---|
| 212 | fid.createDimension('numbers_in_range', 2) |
---|
| 213 | |
---|
| 214 | fid.createVariable('x', 'd', ('number_of_points',)) |
---|
| 215 | fid.createVariable('y', 'd', ('number_of_points',)) |
---|
| 216 | fid.createVariable('elevation', 'd', ('number_of_points',)) |
---|
| 217 | fid.createVariable('elevation_range', 'd', ('numbers_in_range',)) |
---|
| 218 | fid.createVariable('time', 'i', ('number_of_timesteps',)) |
---|
| 219 | fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points')) |
---|
| 220 | fid.createVariable('stage_range', 'd', ('numbers_in_range', )) |
---|
| 221 | fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points')) |
---|
| 222 | fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',)) |
---|
| 223 | fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points')) |
---|
| 224 | fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',)) |
---|
| 225 | |
---|
| 226 | fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
| 227 | fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
| 228 | fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
| 229 | fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036]) |
---|
| 230 | fid.variables['elevation'][:] = elevation |
---|
[7692] | 231 | fid.variables['time'][:] = time#[j*60 : j*60 + 60] |
---|
[7678] | 232 | |
---|
| 233 | s = fid.variables['stage'] |
---|
| 234 | xm = fid.variables['xmomentum'] |
---|
| 235 | ym = fid.variables['ymomentum'] |
---|
| 236 | |
---|
[7794] | 237 | #jj = j*number_of_points*(number_of_timesteps/2) |
---|
| 238 | |
---|
[7678] | 239 | for i in xrange(number_of_timesteps): |
---|
[7692] | 240 | ii = i*number_of_points# + jj |
---|
[7678] | 241 | s[i] = stage[ii : ii + number_of_points] |
---|
| 242 | xm[i] = xmomentum[ii : ii + number_of_points] |
---|
| 243 | ym[i] = ymomentum[ii : ii + number_of_points] |
---|
[7692] | 244 | |
---|
[7678] | 245 | origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings - it works for long and lat |
---|
| 246 | geo_ref = write_NetCDF_georeference(origin, fid) |
---|
| 247 | |
---|
| 248 | fid.variables['x'][:] = x - geo_ref.get_xllcorner() |
---|
| 249 | fid.variables['y'][:] = y - geo_ref.get_yllcorner() |
---|
| 250 | |
---|
[7794] | 251 | fid.close() |
---|
| 252 | |
---|