diff --git a/oggm/core/gis.py b/oggm/core/gis.py index ea9ad9c35..3d81a475c 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -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: @@ -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': @@ -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, @@ -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) @@ -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), @@ -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') @@ -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' @@ -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', ), @@ -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', ), @@ -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' diff --git a/oggm/graphics.py b/oggm/graphics.py index 1241aa3f0..57a4f023a 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..13b9f0ef6 100644 --- a/oggm/sandbox/distribute_2d.py +++ b/oggm/sandbox/distribute_2d.py @@ -9,7 +9,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 from oggm.utils import ncDataset, entity_task import matplotlib.pyplot as plt @@ -416,3 +416,53 @@ def distribute_thickness_from_simulation(gdir, return ds, out_df return ds + + +def combine_distributed_thickness(gdirs, output_folder, 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) + + #create gridded_data.nc - maybe use this further down, where we can already + # add the Topography to the gridded_data + with GriddedNcdfFile(grid=combined_grid, fpath=output_folder) as nc: + # need to check which other data is necessary for the plotting of the + # distributed_thickness? or which other + # data should be included? glacier outlines etc.? + pass + + # retrieve and project a topography/DEM for the defined region + get_dem_for_grid(combined_grid, output_folder, source=source) + + # add 'write_dem_to_gridded_data' function and add dem to the gridded data + + + # add individual add individual distributed thicknesses to gridded data + # for gdir in gdirs: + # with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: + # ds = ds.load() + # + # r_data = combined_grid.map_gridded_data(ds[f'simulation_distributed_thickness_{ssp}_projection'], + # grid=gdir.grid) + # total_data += r_data.filled(0) + # total_data = np.where(total_data == 0, np.NaN, total_data) + + diff --git a/oggm/utils/_funcs.py b/oggm/utils/_funcs.py index 0470a91b2..e3b13c729 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.