[7678] | 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() |
---|