diff --git a/oggm/core/gis.py b/oggm/core/gis.py index 0de25f321..cc91d2776 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -6,6 +6,7 @@ import os import logging import warnings +import json from packaging.version import Version from functools import partial @@ -301,82 +302,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 - """ - - utm_proj, nx, ny, ulx, uly, dx = glacier_grid_params(gdir) + 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. - # 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), - gdir=gdir, - 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: @@ -395,19 +393,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': @@ -418,12 +421,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, @@ -433,7 +433,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) @@ -442,6 +442,119 @@ def _get_nodata(rio_ds): for dem_ds in dem_dss: dem_ds.close() + +def get_dem_for_grid(grid, fpath, source=None, gdir=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]] + if gdir is not None: + rgi_id = gdir.rgi_id + else: + rgi_id = None + + 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), + gdir=gdir, + 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, gdir=gdir) + # 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), @@ -515,19 +628,34 @@ def rasterio_to_gdir(gdir, input_file, output_file_name, dst.write(dest, indexes=i) -def read_geotiff_dem(gdir): - """Reads (and masks out) the DEM out of the gdir's geotiff file. +def read_geotiff_dem(gdir=None, fpath=None): + """Reads (and masks out) the DEM out of the gdir's geotiff file or a geotiff file given by the fpath variable. Parameters ---------- gdir : :py:class:`oggm.GlacierDirectory` the glacier directory + fpath : 'str' + path to a gdir Returns ------- 2D np.float32 array """ - with rasterio.open(gdir.get_filepath('dem'), 'r', driver='GTiff') as ds: + if gdir is not None: + dem_path = gdir.get_filepath('dem') + else: + if fpath is not None: + if fpath[-4:] != '.tif': + # we add the default filename + fpath = os.path.join(fpath, 'dem.tif') + dem_path = fpath + else: + raise InvalidParamsError('If you do not provide a gdir you must' + f'define a fpath! Given fpath={fpath}.') + + + with rasterio.open(dem_path, 'r', driver='GTiff') as ds: topo = ds.read(1).astype(rasterio.float32) topo[topo <= -999.] = np.NaN topo[ds.read_masks(1) == 0] = np.NaN @@ -540,9 +668,38 @@ 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 extent 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 + '.nc') + self.grid = grid + if reset and os.path.exists(self.fpath): os.remove(self.fpath) @@ -586,7 +743,7 @@ def __exit__(self, exc_type, exc_value, exc_traceback): @entity_task(log, writes=['gridded_data']) -def process_dem(gdir): +def process_dem(gdir=None, grid=None, fpath=None): """Reads the DEM from the tiff, attempts to fill voids and apply smooth. The data is then written to `gridded_data.nc`. @@ -596,13 +753,39 @@ def process_dem(gdir): gdir : :py:class:`oggm.GlacierDirectory` where to write the data """ - - # open srtm tif-file: - dem = read_geotiff_dem(gdir) - - # Grid - nx = gdir.grid.nx - ny = gdir.grid.ny + if gdir is not None: + # open srtm tif-file: + dem = read_geotiff_dem(gdir) + # Grid + dem_grid = gdir.grid + grid_name = gdir.rgi_id + 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}.') + dem = read_geotiff_dem(fpath=fpath) + grid_name = 'custom_grid' + dem_grid = grid + diagnostics_dict = dict() + # Grid parameters + nx = dem_grid.nx + ny = dem_grid.ny + dx = dem_grid.dx + xx, yy = dem_grid.ij_coordinates + + def _log_and_diagnostics(pnan, log_string='extra'): + if log_string != 'extra': + perc_string = 'dem_invalid_perc' + else: + perc_string = 'dem_extrapol_perc' + log.info(grid_name + f': DEM needed {log_string}polation.') + if gdir is not None: + gdir.add_to_diagnostics(f'dem_needed_{log_string}polation', True) + gdir.add_to_diagnostics(perc_string, len(pnan[0]) / (nx * ny)) + else: + diagnostics_dict[f'dem_needed_{log_string}polation'] = True + diagnostics_dict[perc_string] = len(pnan[0]) / (nx * ny) # Correct the DEM valid_mask = np.isfinite(dem) @@ -612,8 +795,7 @@ def process_dem(gdir): if np.any(~valid_mask): # We interpolate if np.sum(~valid_mask) > (0.25 * nx * ny): - log.info('({}) more than 25% NaNs in DEM'.format(gdir.rgi_id)) - xx, yy = gdir.grid.ij_coordinates + log.info('({}) more than 25% NaNs in DEM'.format(grid_name)) pnan = np.nonzero(~valid_mask) pok = np.nonzero(valid_mask) points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T @@ -623,15 +805,13 @@ def process_dem(gdir): method='linear') except ValueError: raise InvalidDEMError('DEM interpolation not possible.') - log.info(gdir.rgi_id + ': DEM needed interpolation.') - gdir.add_to_diagnostics('dem_needed_interpolation', True) - gdir.add_to_diagnostics('dem_invalid_perc', len(pnan[0]) / (nx * ny)) + + _log_and_diagnostics(pnan, log_string='inter') isfinite = np.isfinite(dem) if np.any(~isfinite): # interpolation will still leave NaNs in DEM: # extrapolate with NN if needed (e.g. coastal areas) - xx, yy = gdir.grid.ij_coordinates pnan = np.nonzero(~isfinite) pok = np.nonzero(isfinite) points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T @@ -641,13 +821,11 @@ def process_dem(gdir): method='nearest') except ValueError: raise InvalidDEMError('DEM extrapolation not possible.') - log.info(gdir.rgi_id + ': DEM needed extrapolation.') - gdir.add_to_diagnostics('dem_needed_extrapolation', True) - gdir.add_to_diagnostics('dem_extrapol_perc', len(pnan[0]) / (nx * ny)) + _log_and_diagnostics(pnan, log_string='extra') if np.min(dem) == np.max(dem): raise InvalidDEMError('({}) min equal max in the DEM.' - .format(gdir.rgi_id)) + .format(grid_name)) # Clip topography to 0 m a.s.l. if cfg.PARAMS['clip_dem_to_zero']: @@ -655,7 +833,7 @@ def process_dem(gdir): # Smooth DEM? if cfg.PARAMS['smooth_window'] > 0.: - gsize = np.rint(cfg.PARAMS['smooth_window'] / gdir.grid.dx) + gsize = np.rint(cfg.PARAMS['smooth_window'] / dx) smoothed_dem = gaussian_blur(dem, int(gsize)) else: smoothed_dem = dem.copy() @@ -664,8 +842,12 @@ def process_dem(gdir): if cfg.PARAMS['clip_dem_to_zero']: utils.clip_min(smoothed_dem, 0, out=smoothed_dem) + if gdir is None: + with open(os.path.join(fpath, 'dem_diagnostics.json'), 'w') as f: + json.dump(diagnostics_dict, f) + # Write to file - with GriddedNcdfFile(gdir, reset=True) as nc: + with GriddedNcdfFile(gdir=gdir, grid=dem_grid, fpath=fpath, reset=True) as nc: v = nc.createVariable('topo', 'f4', ('y', 'x',), zlib=True) v.units = 'm' @@ -774,7 +956,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', ), diff --git a/oggm/graphics.py b/oggm/graphics.py index 25f87b595..af4bf5023 100644 --- a/oggm/graphics.py +++ b/oggm/graphics.py @@ -76,79 +76,6 @@ def surf_to_nan(surf_h, thick): return surf_h -def combine_grids(gdirs): - """ Combines individual grids of different glacier directories to show - multiple glaciers in the same plot. The resulting grid extent includes - all individual grids completely. - - Parameters - ---------- - gdirs : [], required - A list of GlacierDirectories. - - Returns - ------- - salem.gis.Grid - """ - - new_grid = { - 'proj': None, - 'nxny': None, - 'dxdy': None, - 'x0y0': None, - 'pixel_ref': None - } - - left_use = None - right_use = None - bottom_use = None - top_use = None - dx_use = None - dy_use = None - - for gdir in gdirs: - # use the first gdir to define some values - if new_grid['proj'] is None: - new_grid['proj'] = gdir.grid.proj - if new_grid['pixel_ref'] is None: - new_grid['pixel_ref'] = gdir.grid.pixel_ref - - # find largest extend including all grids completely - (left, right, bottom, top) = gdir.grid.extent_in_crs(new_grid['proj']) - if (left_use is None) or (left_use > left): - left_use = left - if right_use is None or right_use < right: - right_use = right - if bottom_use is None or bottom_use > bottom: - bottom_use = bottom - if top_use is None or top_use < top: - top_use = top - - # find smallest dx and dy for the estimation of nx and ny - dx = gdir.grid.dx - dy = gdir.grid.dy - if dx_use is None or dx_use > dx: - dx_use = dx - # dy could be negative - if dy_use is None or abs(dy_use) > abs(dy): - dy_use = dy - - # calculate nx and ny, the final extend could be one grid point larger or - # smaller due to round() - nx_use = round((right_use - left_use) / dx_use) - ny_use = round((top_use - bottom_use) / abs(dy_use)) - - # finally define the last values of the new grid - if np.sign(dy_use) < 0: - new_grid['x0y0'] = (left_use, top_use) - else: - new_grid['x0y0'] = (left_use, bottom_use) - new_grid['nxny'] = (nx_use, ny_use) - new_grid['dxdy'] = (dx_use, dy_use) - - return salem.gis.Grid.from_dict(new_grid) - - def _plot_map(plotfunc): """ Decorator for common salem.Map plotting logic @@ -213,7 +140,7 @@ def newplotfunc(gdirs, ax=None, smap=None, add_colorbar=True, title=None, if smap is None: if extend_plot_limit: - grid_combined = combine_grids(gdirs) + grid_combined = utils.combine_grids(gdirs) mp = salem.Map(grid_combined, countries=False, nx=grid_combined.nx) else: diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py index d8fdcdb76..3e47ef9ea 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -1,6 +1,7 @@ import logging import warnings +import os import oggm import oggm.cfg as cfg from oggm import utils @@ -9,7 +10,7 @@ import xarray as xr from scipy import ndimage from scipy.stats import mstats -from oggm.core.gis import gaussian_blur +from oggm.core.gis import gaussian_blur, get_dem_for_grid, GriddedNcdfFile, process_dem from oggm.utils import ncDataset, entity_task import matplotlib.pyplot as plt @@ -416,3 +417,107 @@ def distribute_thickness_from_simulation(gdir, return ds, out_df return ds + + +def combine_distributed_thickness(gdirs, output_folder, suffix='', source=None): + """ + This function takes a list of glacier directories that have a + distributed_thickness. + + Parameters + ---------- + gdirs : list of :py:class:`oggm.GlacierDirectory` objects + The glacier directories which should be combined + output_folder : str + Folder where the intermediate files and the final combined gridded data + should be stored + source : str + DEM source + + Returns + ------- + + """ + # create a combined salem.Grid object, which serves as canvas/boundaries of + # the combined glacier region + combined_grid = utils.combine_grids(gdirs) + + # check if output_folder exists + if not os.path.exists(output_folder): + os.makedirs(output_folder) + + # retrieve and project a topography/DEM for the defined region + get_dem_for_grid(combined_grid, output_folder, source=source) + # processes dem(smoothing etc.) and creates gridded_data and adds DEM data + process_dem.unwrapped(gdir=None, grid=combined_grid, fpath=output_folder) + + + def _combine_variable_on_common_grid(var_name, source_file='gridded_data', time_dim=False, suffix=''): + combined_data = None + time_coord = None + for gdir in gdirs: + with xr.open_dataset(gdir.get_filepath(source_file, filesuffix=suffix)) as ds: + ds = ds.load() + + # assuming that the same simulation ran the same time for each glacier + # => taking the simulation duration from the first gdir + if combined_data is None: + if time_dim: + if time_coord is None: + time_coord = ds.time + else: + if time_coord == ds.time: + continue + else: + raise ValueError('The simulation times of the simulations that you try to combine are ' + 'different from each other.\n' + f'The Glacier with the rgi_id {gdir.rgi_id}\ has {len(ds.time)} timesteps,' + f'not matching with {gdirs[0].rgi_id}\'s {len(time_coord)} timesteps.') + combined_data = np.zeros((len(time_coord), combined_grid.ny, combined_grid.nx)) + else: + combined_data = np.zeros((combined_grid.ny, combined_grid.nx)) + r_data = combined_grid.map_gridded_data(ds[var_name], + grid=gdir.grid) + combined_data += r_data.filled(0) + combined_data = np.where(combined_data == 0, np.NaN, combined_data) + + return combined_data, time_coord + + combined_distributed_thickness, time_coord = _combine_variable_on_common_grid('distributed_thickness', + source_file='gridded_simulation', time_dim=True, + suffix=suffix) + combined_glacier_mask, _ = _combine_variable_on_common_grid('glacier_mask') + # combined_topo_smoothed = _combine_variable_on_common_grid('topo_smoothed') + combined_bedrock, _ = _combine_variable_on_common_grid('bedrock', source_file='gridded_simulation', suffix=suffix) + combined_glacier_ext, _ = _combine_variable_on_common_grid('glacier_ext') + + + with utils.ncDataset(os.path.join(output_folder, 'gridded_data.nc'), 'a') as gridded_data: + gridded_data.createDimension('time', len(time_coord)) + v = gridded_data.createVariable('time', 'f4', ('time',), zlib=True) + v.units = 'years' + v.long_name = 'time coordinate of the simulation' + v.standard_name = 'time' + v[:] = time_coord.values + # add combined distributed thickness + + v = gridded_data.createVariable('combined_distributed_thickness', 'f4', ('time', 'y', 'x',), zlib=True) + v.units = 'm' + v.long_name = 'combined distributed thickness' + v[:] = combined_distributed_thickness + + v = gridded_data.createVariable('combined_glacier_mask', 'f4', ('y', 'x',), zlib=True) + v.units = 'm' + v.long_name = 'combined glacier mask' + v[:] = combined_glacier_mask + + v = gridded_data.createVariable('combined_bedrock', 'f4', ('y', 'x',), zlib=True) + v.units = 'm' + v.long_name = 'combined bedrock' + v[:] = combined_bedrock + + v = gridded_data.createVariable('combined_glacier_ext', 'f4', ('y', 'x',), zlib=True) + v.units = 'm' + v.long_name = 'combined glacier extent' + v[:] = combined_glacier_ext + diff --git a/oggm/utils/_funcs.py b/oggm/utils/_funcs.py index 644fe179a..cffd11d97 100644 --- a/oggm/utils/_funcs.py +++ b/oggm/utils/_funcs.py @@ -24,6 +24,7 @@ import shapely.geometry as shpg from shapely.ops import linemerge from shapely.validation import make_valid +from salem.gis import Grid # Optional libs try: @@ -539,6 +540,79 @@ def recursive_valid_polygons(geoms, crs): return new_geoms +def combine_grids(gdirs): + """ Combines individual grids of different glacier directories to show + multiple glaciers in the same plot. The resulting grid extent includes + all individual grids completely. + + Parameters + ---------- + gdirs : [], required + A list of GlacierDirectories. + + Returns + ------- + salem.gis.Grid + """ + + new_grid = { + 'proj': None, + 'nxny': None, + 'dxdy': None, + 'x0y0': None, + 'pixel_ref': None + } + + left_use = None + right_use = None + bottom_use = None + top_use = None + dx_use = None + dy_use = None + + for gdir in gdirs: + # use the first gdir to define some values + if new_grid['proj'] is None: + new_grid['proj'] = gdir.grid.proj + if new_grid['pixel_ref'] is None: + new_grid['pixel_ref'] = gdir.grid.pixel_ref + + # find largest extend including all grids completely + (left, right, bottom, top) = gdir.grid.extent_in_crs(new_grid['proj']) + if (left_use is None) or (left_use > left): + left_use = left + if right_use is None or right_use < right: + right_use = right + if bottom_use is None or bottom_use > bottom: + bottom_use = bottom + if top_use is None or top_use < top: + top_use = top + + # find smallest dx and dy for the estimation of nx and ny + dx = gdir.grid.dx + dy = gdir.grid.dy + if dx_use is None or dx_use > dx: + dx_use = dx + # dy could be negative + if dy_use is None or abs(dy_use) > abs(dy): + dy_use = dy + + # calculate nx and ny, the final extend could be one grid point larger or + # smaller due to round() + nx_use = round((right_use - left_use) / dx_use) + ny_use = round((top_use - bottom_use) / abs(dy_use)) + + # finally define the last values of the new grid + if np.sign(dy_use) < 0: + new_grid['x0y0'] = (left_use, top_use) + else: + new_grid['x0y0'] = (left_use, bottom_use) + new_grid['nxny'] = (nx_use, ny_use) + new_grid['dxdy'] = (dx_use, dy_use) + + return Grid.from_dict(new_grid) + + def multipolygon_to_polygon(geometry, gdir=None): """Sometimes an RGI geometry is a multipolygon: this should not happen. diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py index 88f49366a..df75540d3 100644 --- a/oggm/utils/_workflow.py +++ b/oggm/utils/_workflow.py @@ -515,6 +515,8 @@ def _entity_task(gdir, *, reset=None, print_log=True, return out _entity_task.__dict__['is_entity_task'] = True + # adds the possibility to use a function, decorated as entity_task, without its decoration. + _entity_task.unwrapped = task_func return _entity_task