Skip to content

Small adaptations to your code #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pat-schmitt wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Small adaptations to your code #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
a4a1e01
moved combine_grids function from graphics to utils
afisc Nov 27, 2023
816cc25
introduced 'combine_distributed_thickness' functionality - not finish…
afisc Nov 27, 2023
2296a14
split up and redesigned the 'define_glacier_region' to also use parts…
afisc Nov 27, 2023
874b2d0
bugfix
afisc Nov 27, 2023
f764bab
bugfix
afisc Nov 28, 2023
628ef97
added comments to the combine_distributed_thickness function
afisc Nov 28, 2023
83781d4
cleaning and refractoring
pat-schmitt Nov 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
291 changes: 213 additions & 78 deletions oggm/core/gis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -298,82 +298,79 @@ def glacier_grid_params(gdir):
return utm_proj, nx, ny, ulx, uly, dx


@entity_task(log, writes=['glacier_grid', 'dem', 'outlines'])
def define_glacier_region(gdir, entity=None, source=None):
"""Very first task after initialization: define the glacier's local grid.

Defines the local projection (Transverse Mercator or UTM depending on
user choice), centered on the glacier.
There is some options to set the resolution of the local grid.
It can be adapted depending on the size of the glacier.
See ``params.cfg`` for setting these options.

Default values of the adapted mode lead to a resolution of 50 m for
Hintereisferner, which is approx. 8 km2 large.

After defining the grid, the topography and the outlines of the glacier
are transformed into the local projection.
def check_dem_source(source, extent_ll, rgi_id=None):
"""
This function can check for multiple DEM sources and is in charge of the
error handling if a requested source is not available for the given
glacier/extent

Parameters
----------
gdir : :py:class:`oggm.GlacierDirectory`
where to write the data
entity : geopandas.GeoSeries
the glacier geometry to process - DEPRECATED. It is now ignored
source : str or list of str, optional
If you want to force the use of a certain DEM source. Available are:
- 'USER' : file set in cfg.PATHS['dem_file']
- 'SRTM' : http://srtm.csi.cgiar.org/
- 'GIMP' : https://bpcrc.osu.edu/gdg/data/gimpdem
- 'RAMP' : https://nsidc.org/data/nsidc-0082/versions/2/documentation
- 'REMA' : https://www.pgc.umn.edu/data/rema/
- 'DEM3' : http://viewfinderpanoramas.org/
- 'ASTER' : https://asterweb.jpl.nasa.gov/gdem.asp
- 'TANDEM' : https://geoservice.dlr.de/web/dataguide/tdm90/
- 'ARCTICDEM' : https://www.pgc.umn.edu/data/arcticdem/
- 'AW3D30' : https://www.eorc.jaxa.jp
- 'MAPZEN' : https://registry.opendata.aws/terrain-tiles/
- 'ALASKA' : https://www.the-cryosphere.net/8/503/2014/
- 'COPDEM30' : Copernicus DEM GLO30 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
- 'COPDEM90' : Copernicus DEM GLO90 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
- 'NASADEM': https://doi.org/10.5069/G93T9FD9
"""
source : str or list of str
If you want to force the use of a certain DEM source. If list is
provided they are checked in order. For a list of available options
check docstring of oggm.core.gis.define_glacier_region.
extent_ll : list
The longitude and latitude extend which should be checked for. Should
be provided as extent_ll = [[minlon, maxlon], [minlat, maxlat]].
rgi_id : str
The RGI-ID of the glacier, only used for a descriptive error message in
case.

utm_proj, nx, ny, ulx, uly, dx = glacier_grid_params(gdir)

# Back to lon, lat for DEM download/preparation
tmp_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), x0y0=(ulx, uly),
dxdy=(dx, -dx), pixel_ref='corner')
minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84)
Returns
-------
String of the working DEM source.
"""

# Open DEM
# We test DEM availability for glacier only (maps can grow big)
if isinstance(source, list): # when multiple sources are provided, try them sequentially
# when multiple sources are provided, try them sequentially
if isinstance(source, list):
for src in source:
source_exists = is_dem_source_available(src, *gdir.extent_ll)
source_exists = is_dem_source_available(src, *extent_ll)
if source_exists:
source = src # pick the first source which exists
break
else:
source_exists = is_dem_source_available(source, *gdir.extent_ll)

source_exists = is_dem_source_available(source, *extent_ll)
if not source_exists:
raise InvalidWorkflowError(f'Source: {source} is not available for '
f'glacier {gdir.rgi_id} with border '
f"{cfg.PARAMS['border']}")
dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
rgi_id=gdir.rgi_id,
dx_meter=dx,
source=source)
log.debug('(%s) DEM source: %s', gdir.rgi_id, dem_source)
log.debug('(%s) N DEM Files: %s', gdir.rgi_id, len(dem_list))
if rgi_id is None:
extent_string = (f"the grid extent of longitudes {extent_ll[0]} "
f"and latitudes {extent_ll[1]}")
else:
extent_string = (f"the glacier {rgi_id} with border "
f"{cfg.PARAMS['border']}")
raise InvalidWorkflowError(f"Source: {source} is not available for "
f"{extent_string}")
return source


def reproject_dem(dem_list, dem_source, dst_grid_prop, output_path):
"""
This function reprojects a provided DEM to the destination grid and saves
the result to disk.

Parameters
----------
dem_list : list of str
list with path(s) to the DEM file(s)
dem_source : str
DEM source string
dst_grid_prop : dict
Holds necessary grid properties. Must contain 'utm_proj', 'dx', 'ulx',
'uly', 'nx' and 'ny.
output_path : str
Filepath where to store the reprojected DEM.

Returns
-------
None
"""

# Decide how to tag nodata
def _get_nodata(rio_ds):
nodata = rio_ds[0].meta.get('nodata', None)
if nodata is None:
# badly tagged geotiffs, let's do it ourselves
nodata = -32767 if source == 'TANDEM' else -9999
nodata = -32767 if dem_source == 'TANDEM' else -9999
return nodata

# A glacier area can cover more than one tile:
Expand All @@ -392,19 +389,24 @@ def _get_nodata(rio_ds):

# Use Grid properties to create a transform (see rasterio cookbook)
dst_transform = rasterio.transform.from_origin(
ulx, uly, dx, dx # sign change (2nd dx) is done by rasterio.transform
dst_grid_prop['ulx'], dst_grid_prop['uly'], dst_grid_prop['dx'],
dst_grid_prop['dx']
# sign change (2nd dx) is done by rasterio.transform
)

# Set up profile for writing output
profile = dem_dss[0].profile
profile.update({
'crs': utm_proj.srs,
'crs': dst_grid_prop['utm_proj'].srs,
'transform': dst_transform,
'nodata': nodata,
'width': nx,
'height': ny,
'width': dst_grid_prop['nx'],
'height': dst_grid_prop['ny'],
'driver': 'GTiff'
})
profile.pop('blockxsize', None)
profile.pop('blockysize', None)
profile.pop('compress', None)

# Could be extended so that the cfg file takes all Resampling.* methods
if cfg.PARAMS['topo_interp'] == 'bilinear':
Expand All @@ -415,12 +417,9 @@ def _get_nodata(rio_ds):
raise InvalidParamsError('{} interpolation not understood'
.format(cfg.PARAMS['topo_interp']))

dem_reproj = gdir.get_filepath('dem')
profile.pop('blockxsize', None)
profile.pop('blockysize', None)
profile.pop('compress', None)
with rasterio.open(dem_reproj, 'w', **profile) as dest:
dst_array = np.empty((ny, nx), dtype=dem_dss[0].dtypes[0])
with rasterio.open(output_path, 'w', **profile) as dest:
dst_array = np.empty((dst_grid_prop['ny'], dst_grid_prop['nx']),
dtype=dem_dss[0].dtypes[0])
reproject(
# Source parameters
source=dem_data,
Expand All @@ -430,7 +429,7 @@ def _get_nodata(rio_ds):
# Destination parameters
destination=dst_array,
dst_transform=dst_transform,
dst_crs=utm_proj.srs,
dst_crs=dst_grid_prop['utm_proj'].srs,
dst_nodata=nodata,
# Configuration
resampling=resampling)
Expand All @@ -439,6 +438,114 @@ def _get_nodata(rio_ds):
for dem_ds in dem_dss:
dem_ds.close()


def get_dem_for_grid(grid, fpath, source=None, rgi_id=None):
"""
Fetch a DEM from source, reproject it to the extent defined by grid and
saves it to disk.

Parameters
----------
grid : :py:class:`salem.gis.Grid`
Grid which defines the extent and projection of the final DEM
fpath : str
The output filepath for the final DEM.
source : str or list of str
If you want to force the use of a certain DEM source. If list is
provided they are checked in order. For a list of available options
check docstring of oggm.core.gis.define_glacier_region.
rgi_id : str
The RGI-ID of the glacier.

Returns
-------
tuple: (list with path(s) to the DEM file(s), data source str)
"""
minlon, maxlon, minlat, maxlat = grid.extent_in_crs(crs=salem.wgs84)
extent_ll = [[minlon, maxlon], [minlat, maxlat]]
grid_prop = {
'utm_proj': grid.proj,
'dx': grid.dx,
'ulx': grid.x0,
'uly': grid.y0,
'nx': grid.nx,
'ny': grid.ny
}

source = check_dem_source(source, extent_ll, rgi_id=rgi_id)

dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
rgi_id=rgi_id,
dx_meter=grid_prop['dx'],
source=source)

if rgi_id is not None:
log.debug('(%s) DEM source: %s', rgi_id, dem_source)
log.debug('(%s) N DEM Files: %s', rgi_id, len(dem_list))

# further checks if given fpath exists?
if fpath[-4:] != '.tif':
# we add the default filename
fpath = os.path.join(fpath, 'dem.tif')

reproject_dem(dem_list=dem_list, dem_source=dem_source,
dst_grid_prop=grid_prop,
output_path=fpath)

return dem_list, dem_source


@entity_task(log, writes=['glacier_grid', 'dem', 'outlines'])
def define_glacier_region(gdir, entity=None, source=None):
"""Very first task after initialization: define the glacier's local grid.

Defines the local projection (Transverse Mercator or UTM depending on
user choice), centered on the glacier.
There is some options to set the resolution of the local grid.
It can be adapted depending on the size of the glacier.
See ``params.cfg`` for setting these options.

Default values of the adapted mode lead to a resolution of 50 m for
Hintereisferner, which is approx. 8 km2 large.

After defining the grid, the topography and the outlines of the glacier
are transformed into the local projection.

Parameters
----------
gdir : :py:class:`oggm.GlacierDirectory`
where to write the data
entity : geopandas.GeoSeries
the glacier geometry to process - DEPRECATED. It is now ignored
source : str or list of str, optional
If you want to force the use of a certain DEM source. Available are:
- 'USER' : file set in cfg.PATHS['dem_file']
- 'SRTM' : http://srtm.csi.cgiar.org/
- 'GIMP' : https://bpcrc.osu.edu/gdg/data/gimpdem
- 'RAMP' : https://nsidc.org/data/nsidc-0082/versions/2/documentation
- 'REMA' : https://www.pgc.umn.edu/data/rema/
- 'DEM3' : http://viewfinderpanoramas.org/
- 'ASTER' : https://asterweb.jpl.nasa.gov/gdem.asp
- 'TANDEM' : https://geoservice.dlr.de/web/dataguide/tdm90/
- 'ARCTICDEM' : https://www.pgc.umn.edu/data/arcticdem/
- 'AW3D30' : https://www.eorc.jaxa.jp
- 'MAPZEN' : https://registry.opendata.aws/terrain-tiles/
- 'ALASKA' : https://www.the-cryosphere.net/8/503/2014/
- 'COPDEM30' : Copernicus DEM GLO30 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
- 'COPDEM90' : Copernicus DEM GLO90 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
- 'NASADEM': https://doi.org/10.5069/G93T9FD9
"""

utm_proj, nx, ny, ulx, uly, dx = glacier_grid_params(gdir)

# Back to lon, lat for DEM download/preparation
tmp_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), x0y0=(ulx, uly),
dxdy=(dx, -dx), pixel_ref='corner')

dem_list, dem_source = get_dem_for_grid(grid=tmp_grid,
fpath=gdir.get_filepath('dem'),
source=source, rgi_id=gdir.rgi_id)

# Glacier grid
x0y0 = (ulx+dx/2, uly-dx/2) # To pixel center coordinates
glacier_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), dxdy=(dx, -dx),
Expand Down Expand Up @@ -537,14 +644,42 @@ class GriddedNcdfFile(object):
The other variables have to be created and filled by the calling
routine.
"""
def __init__(self, gdir, basename='gridded_data', reset=False):
self.fpath = gdir.get_filepath(basename)
self.grid = gdir.grid
def __init__(self, gdir=None, grid=None, fpath=None,
basename='gridded_data', reset=False):
"""
Parameters
----------
gdir : :py:class:`oggm.GlacierDirectory`
The glacier directory. If provided it defines the filepath and the
grid used for the netcdf file. It overrules grid and fpath if all
are provided.
grid : :py:class:`salem.gis.Grid`
Grid which defines the extend of the netcdf file. Needed if gdir is
not provided.
fpath : str
The output filepath for the netcdf file. Needed if gdir is not
provided.
basename : str
The filename of the resulting netcdf file
reset : bool
If True, a potentially existing file will be deleted.
"""

if gdir is not None:
self.fpath = gdir.get_filepath(basename)
self.grid = gdir.grid
else:
if grid is None or fpath is None:
raise InvalidParamsError('If you do not provide a gdir you must'
'define grid and fpath! Given grid='
f'{grid} and fpath={fpath}.')
self.fpath = os.path.join(fpath, basename)
self.grid = grid

if reset and os.path.exists(self.fpath):
os.remove(self.fpath)

def __enter__(self):

if os.path.exists(self.fpath):
# Already there - just append
self.nc = ncDataset(self.fpath, 'a', format='NETCDF4')
Expand Down Expand Up @@ -662,7 +797,7 @@ def process_dem(gdir):
utils.clip_min(smoothed_dem, 0, out=smoothed_dem)

# Write to file
with GriddedNcdfFile(gdir, reset=True) as nc:
with GriddedNcdfFile(gdir=gdir, reset=True) as nc:

v = nc.createVariable('topo', 'f4', ('y', 'x',), zlib=True)
v.units = 'm'
Expand Down Expand Up @@ -771,7 +906,7 @@ def proj(x, y):
gdir.write_pickle(geometries, 'geometries')

# write out the grids in the netcdf file
with GriddedNcdfFile(gdir) as nc:
with GriddedNcdfFile(gdir=gdir) as nc:

if 'glacier_mask' not in nc.variables:
v = nc.createVariable('glacier_mask', 'i1', ('y', 'x', ),
Expand Down Expand Up @@ -905,7 +1040,7 @@ def simple_glacier_masks(gdir):
.format(gdir.rgi_id))

# write out the grids in the netcdf file
with GriddedNcdfFile(gdir) as nc:
with GriddedNcdfFile(gdir=gdir) as nc:

if 'glacier_mask' not in nc.variables:
v = nc.createVariable('glacier_mask', 'i1', ('y', 'x', ),
Expand Down Expand Up @@ -1619,7 +1754,7 @@ def project(x, y):
.format(gdir.rgi_id))

# write out the grids in the netcdf file
with GriddedNcdfFile(gdir, reset=True) as nc:
with GriddedNcdfFile(gdir=gdir, reset=True) as nc:

v = nc.createVariable('topo', 'f4', ('y', 'x', ), zlib=True)
v.units = 'm'
Expand Down
Loading