Often we end up with several versions of a large topo file that have been cropped/coarsened in different ways for specific regions, or that might have a simple change like shifting x in header by 360 to change from longitude E to W
Rather than doing this in Python and saving multiple versions, in many cases it would be preferable to reuse a single DEM that covers a large area (e.g. etopo or GEBCO data that covers the entire Pacific) and read the entire file into GeoClaw (preferably the netCDF version) and then crop, coarsen, and shift as desired in the Fortran code.
To facilitate this, I've been talking with @mandli and we propose adding more flexibility to the
rundata.topo_data.topofiles
list that is set in setrun.py. Currently this is a list of 2-element lists (or tuples) of the form [topo_type, topo_path].
I suggest allowing the topofiles list to contain either such lists (for backward compatibility) or dictionaries (or a mix of the two). Default values indicated below would be used for any missing dictionary keys.
A dictionary topofile would contain key-value pairs:
- 'topo_type' : as before
- 'topo_path' : as before
- 'extent' : [x1, x2, y1, y2] if cropping is desired, or 'all' (default 'all')
- 'coarsen' : factor to coarsen by (default 1)
- 'buffer' : buffering of extent (as in topotools.Topography.crop()) (default 0)
- 'align' : alignment, if coarsening (as in topotools.Topography.crop()) (default None)
- 'x_shift' : shift to apply to longitude x, e.g. -360. (default 0)
- 'z_shift' : shift to apply to Z, e.g. to adjust vertical datum (default 0)
- 'negate_z' : if True, negate Z before shifting (default False)
After reading the file into temporary storage, it would be cropped, coarsened, and shifted if desired, and then only the result would be stored in a topo array for later use in GeoClaw.
Are there any other operations one might want to apply to the topo DEM?
This will also require changing the information written into topo.data by make data, and that will be read by the Fortran. Currently we write 2 lines for each topofile:
so we could just add more lines for each topofile corresponding to all the keys listed above, filling in the default values as needed. Is there a better way?
By the way, one might wonder if it is necessary to crop topo files, or if it's ok to read in a topofile that covers a much larger region than is needed at that resolution, in order to avoid cropping. For example, if one wants 1/3" topography around Seaside, OR, why not just read in the full 1/9" CUDEM tile(s) that covers this region and use it all in GeoClaw, rather than cropping and coarsening? Unfortunately Seaside lies at exactly latitude 46.0, so one would need two such tiles, each of which covers 0.25 degrees in x and y with an 8100 x 8100 grid and the combined grids cover about 50 km of N-S coastline centered at the site of interest, at 3 m resolution, also extending eastward far into the mountains. Since GeoClaw always constructs piecewise bilinear functions for every topofile provided, and then integrates this function over the computational grid cells to get the cell-averaged topography values B, even in regions where much coarser computational cells are being used, providing such fine topofiles can add greatly to the computational expense. This is particularly true if a kinematic rupture is specified (e.g. a CSZ earthquake that goes on for several minutes), since then the cell-averaged topo values are recomputed every time step as the topography changes with time during the earthquake.
Often we end up with several versions of a large topo file that have been cropped/coarsened in different ways for specific regions, or that might have a simple change like shifting x in header by 360 to change from longitude E to W
Rather than doing this in Python and saving multiple versions, in many cases it would be preferable to reuse a single DEM that covers a large area (e.g. etopo or GEBCO data that covers the entire Pacific) and read the entire file into GeoClaw (preferably the netCDF version) and then crop, coarsen, and shift as desired in the Fortran code.
To facilitate this, I've been talking with @mandli and we propose adding more flexibility to the
rundata.topo_data.topofileslist that is set in
setrun.py. Currently this is a list of 2-element lists (or tuples) of the form[topo_type, topo_path].I suggest allowing the topofiles list to contain either such lists (for backward compatibility) or dictionaries (or a mix of the two). Default values indicated below would be used for any missing dictionary keys.
A dictionary topofile would contain key-value pairs:
After reading the file into temporary storage, it would be cropped, coarsened, and shifted if desired, and then only the result would be stored in a topo array for later use in GeoClaw.
Are there any other operations one might want to apply to the topo DEM?
This will also require changing the information written into
topo.databymake data, and that will be read by the Fortran. Currently we write 2 lines for each topofile:so we could just add more lines for each topofile corresponding to all the keys listed above, filling in the default values as needed. Is there a better way?
By the way, one might wonder if it is necessary to crop topo files, or if it's ok to read in a topofile that covers a much larger region than is needed at that resolution, in order to avoid cropping. For example, if one wants 1/3" topography around Seaside, OR, why not just read in the full 1/9" CUDEM tile(s) that covers this region and use it all in GeoClaw, rather than cropping and coarsening? Unfortunately Seaside lies at exactly latitude 46.0, so one would need two such tiles, each of which covers 0.25 degrees in
xandywith an8100 x 8100grid and the combined grids cover about 50 km of N-S coastline centered at the site of interest, at 3 m resolution, also extending eastward far into the mountains. Since GeoClaw always constructs piecewise bilinear functions for every topofile provided, and then integrates this function over the computational grid cells to get the cell-averaged topography valuesB, even in regions where much coarser computational cells are being used, providing such fine topofiles can add greatly to the computational expense. This is particularly true if a kinematic rupture is specified (e.g. a CSZ earthquake that goes on for several minutes), since then the cell-averaged topo values are recomputed every time step as the topography changes with time during the earthquake.