From 47ef9f3c88b7bc25d8c3130def8b2b5e076bdc21 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 2 Apr 2024 20:17:46 +0100 Subject: [PATCH 1/2] Add utilities for glacier entities in complexes --- oggm/cfg.py | 6 ++- oggm/core/gis.py | 83 ++++++++++++++++++++++++++++++++++++++++ oggm/tasks.py | 1 + oggm/utils/_downloads.py | 26 +++++++++++++ 4 files changed, 114 insertions(+), 2 deletions(-) diff --git a/oggm/cfg.py b/oggm/cfg.py index 75e60d164..b6d8d273d 100644 --- a/oggm/cfg.py +++ b/oggm/cfg.py @@ -9,7 +9,6 @@ import glob from collections import OrderedDict from distutils.util import strtobool -import warnings import numpy as np import pandas as pd @@ -28,7 +27,7 @@ except ImportError: pass -from oggm.exceptions import InvalidParamsError, InvalidWorkflowError +from oggm.exceptions import InvalidParamsError # Local logger log = logging.getLogger(__name__) @@ -299,6 +298,9 @@ def _log_param_change(self, key, value): _doc = "A table containing the Huss&Farinotti 2012 squeezed flowlines." BASENAMES['elevation_band_flowline'] = ('elevation_band_flowline.csv', _doc) +_doc = "The outlines of this glacier complex sub-entities (for RGI7C only!)." +BASENAMES['complex_sub_entities'] = ('complex_sub_entities.shp', _doc) + def set_logging_config(logging_level='INFO'): """Set the global logger parameters. diff --git a/oggm/core/gis.py b/oggm/core/gis.py index 5b6ed5984..4b63f1182 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -1980,3 +1980,86 @@ def reproject_gridded_data_variable_to_grid(gdir, r_data *= factor return r_data + + +@entity_task(log, writes=['gridded_data', 'complex_sub_entities']) +def rgi7g_to_complex(gdir, rgi7g_file=None, rgi7c_to_g_links=None): + """Adds the individual glacier outlines to this glacier complex gdir. + + Also adds a mask to gridded_data.nc with number indicating the + respective subentity index in the shapefile (-1 means no glacier). + + Parameters + ---------- + gdir : :py:class:`oggm.GlacierDirectory` object + the glacier directory to process + rgi7g_file : gpd.GeoDataFrame + the RGI7G file to extract the outlines from (we can read it + from disk but if you give it, this may faster for large number + of glaciers) + rgi7c_to_g_links : gpd.GeoDataFrame + the RGI7G file to extract the outlines from (we can read it + from disk but if you give it, this may faster for large number + of glaciers) + """ + + if not gdir.rgi_version == '70C': + raise InvalidWorkflowError('Needs to be run on a glacier complex!') + + if rgi7g_file is None: + from oggm.utils import get_rgi_region_file + rgi7g_file = get_rgi_region_file(gdir.rgi_region, version='70G') + rgi7g_file = gpd.read_file(rgi7g_file) + if rgi7c_to_g_links is None: + from oggm.utils import get_rgi7c_to_g_links + rgi7c_to_g_links = get_rgi7c_to_g_links(gdir.rgi_region) + + subset = rgi7g_file.loc[rgi7g_file.rgi_id.isin(rgi7c_to_g_links[gdir.rgi_id])] + subset = subset.reset_index(drop=True) + + # Reproject and write + subset = subset.to_crs(gdir.grid.proj.srs) + gdir.write_shapefile(subset, 'complex_sub_entities') + + # OK all good, now the mask + # Load the DEM + with rasterio.open(gdir.get_filepath('dem'), 'r', driver='GTiff') as ds: + data = ds.read(1).astype(rasterio.float32) + profile = ds.profile + + # Initialize the mask with -1s, the same shape as the DEM + mask = np.full(data.shape, -1, dtype=np.int16) + + # Iterate over each polygon, rasterizing it onto the mask + for index, geometry in enumerate(subset.geometry): + # Correct invalid polygons + geometry = geometry.buffer(0) + if not geometry.is_valid: + raise Exception('Invalid geometry.') + + # Compute the glacier mask using rasterio + # Small detour as mask only accepts DataReader objects + profile['dtype'] = 'int16' + profile.pop('nodata', None) + with rasterio.io.MemoryFile() as memfile: + with memfile.open(**profile) as dataset: + dataset.write(data.astype(np.int16)[np.newaxis, ...]) + dem_data = rasterio.open(memfile.name) + masked_dem, _ = riomask(dem_data, [shpg.mapping(geometry)], + filled=False) + glacier_mask = ~masked_dem[0, ...].mask + + # Update the mask: only change -1 to the new index + mask = np.where((mask == -1) & glacier_mask, index, mask) + + grids_file = gdir.get_filepath('gridded_data') + with ncDataset(grids_file, 'a') as nc: + + vn = 'sub_entities' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'i1', ('y', 'x', )) + v.units = '-' + v.long_name = 'Sub-entities glacier mask (number is index)' + v[:] = mask diff --git a/oggm/tasks.py b/oggm/tasks.py index e09a9766f..f9ed3c769 100644 --- a/oggm/tasks.py +++ b/oggm/tasks.py @@ -13,6 +13,7 @@ from oggm.core.gis import gridded_attributes from oggm.core.gis import gridded_mb_attributes from oggm.core.gis import gridded_data_var_to_geotiff +from oggm.core.gis import rgi7g_to_complex from oggm.core.centerlines import compute_centerlines from oggm.core.centerlines import compute_downstream_line from oggm.core.centerlines import compute_downstream_bedshape diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 5b644da0e..6152b2101 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -23,6 +23,7 @@ import ftplib import ssl import tarfile +import json # External libs import pandas as pd @@ -2231,6 +2232,31 @@ def get_rgi_intersects_entities(rgi_ids, version=None): return selection +def get_rgi7c_to_g_links(region, version='70C', reset=False): + """Path to the RGI7 glacier complex to glacier json file. + + This is for version 7C only! + + Parameters + ---------- + region : str + from '01' to '19' + version : str + '70C', defaults to 70C + reset : bool + If True, deletes the RGI directory first and downloads the data + + Returns + ------- + a dictionary containing the links + """ + jfile = get_rgi_region_file(region, version=version, reset=reset) + jfile = jfile.replace('.shp', '-CtoG_links.json') + with open(jfile) as f: + out = json.load(f) + return out + + def is_dem_source_available(source, lon_ex, lat_ex): """Checks if a DEM source is available for your purpose. From 9a5c2294e434b5655f644da2e76165ce61acfd38 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 2 Apr 2024 20:20:13 +0100 Subject: [PATCH 2/2] pep8 --- oggm/core/gis.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/oggm/core/gis.py b/oggm/core/gis.py index 4b63f1182..88d3e4bd0 100644 --- a/oggm/core/gis.py +++ b/oggm/core/gis.py @@ -1986,19 +1986,19 @@ def reproject_gridded_data_variable_to_grid(gdir, def rgi7g_to_complex(gdir, rgi7g_file=None, rgi7c_to_g_links=None): """Adds the individual glacier outlines to this glacier complex gdir. - Also adds a mask to gridded_data.nc with number indicating the + Also adds a mask to gridded_data.nc with number indicating the respective subentity index in the shapefile (-1 means no glacier). Parameters ---------- gdir : :py:class:`oggm.GlacierDirectory` object the glacier directory to process - rgi7g_file : gpd.GeoDataFrame - the RGI7G file to extract the outlines from (we can read it + rgi7g_file : gpd.GeoDataFrame + the RGI7G file to extract the outlines from (we can read it from disk but if you give it, this may faster for large number of glaciers) - rgi7c_to_g_links : gpd.GeoDataFrame - the RGI7G file to extract the outlines from (we can read it + rgi7c_to_g_links : gpd.GeoDataFrame + the RGI7G file to extract the outlines from (we can read it from disk but if you give it, this may faster for large number of glaciers) """