From 0df0eabda4f2193c6a0a90f3bead616fa5ec0b84 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 19 Jul 2022 19:37:27 +0200 Subject: [PATCH 01/11] WIP: add Millan thickness and velocity to the shop --- oggm/shop/millan22.py | 80 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 oggm/shop/millan22.py diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py new file mode 100644 index 000000000..62b62f472 --- /dev/null +++ b/oggm/shop/millan22.py @@ -0,0 +1,80 @@ +import logging + +import numpy as np +import shapely.geometry as shpg + +try: + import salem +except ImportError: + pass + +try: + import geopandas as gpd +except ImportError: + pass + +from oggm import utils +from oggm.exceptions import InvalidWorkflowError + +# Module logger +log = logging.getLogger(__name__) + +default_base_url = 'https://cluster.klima.uni-bremen.de/~oggm/velocities/millan22/' + +_lookup_thickness = None + +def _get_lookup_thickness(): + global _lookup_thickness + if _lookup_thickness is None: + fname = default_base_url + 'millan22_thickness_lookup_shp.zip' + _lookup_thickness = gpd.read_file('zip://' + utils.file_downloader(fname)) + return _lookup_thickness + +@utils.entity_task(log, writes=['gridded_data']) +def add_millan_thickness(gdir): + """Add the Millan 22 thickness data to this glacier directory. + + Parameters + ---------- + gdir ::py:class:`oggm.GlacierDirectory` + the glacier directory to process + """ + + base_url = default_base_url + + gdf = _get_lookup_thickness() + + cp = shpg.Point(gdir.cenlon, gdir.cenlat) + sel = gdf.loc[gdf.contains(cp)] + if len(sel) > 1: + raise NotImplementedError('Multifile Millan not implemented yet') + if len(sel) == 0: + raise InvalidWorkflowError(f'There seems to be no Millan file for this ' + f'glacier: {gdir.rgi_id}') + + url = base_url + sel['path'].iloc[0] + input_file = utils.file_downloader(url) + + dsb = salem.GeoTiff(input_file) + x0, x1, y0, y1 = gdir.grid.extent + dsb.set_subset(corners=((x0, y0), (x1, y1)), crs=gdir.grid.proj, margin=5) + + thick = dsb.get_vardata().astype(np.float64) + thick = gdir.grid.map_gridded_data(thick, dsb.grid, interp='linear').filled(np.nan) + + # We mask zero ice as nodata + thick = np.where(thick == 0, np.NaN, thick) + + # Write + with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc: + + vn = 'millan_ice_thickness' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x', ), zlib=True) + v.units = 'm' + ln = 'Ice thickness from Millan et al. 2022' + v.long_name = ln + v.base_url = base_url + v[:] = thick From 257ca88acef72759a20d43d7c569279fd781fab8 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Wed, 20 Jul 2022 13:39:08 +0200 Subject: [PATCH 02/11] fixes --- oggm/shop/millan22.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index 62b62f472..d23ceb2e2 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -30,6 +30,7 @@ def _get_lookup_thickness(): _lookup_thickness = gpd.read_file('zip://' + utils.file_downloader(fname)) return _lookup_thickness + @utils.entity_task(log, writes=['gridded_data']) def add_millan_thickness(gdir): """Add the Millan 22 thickness data to this glacier directory. @@ -40,10 +41,8 @@ def add_millan_thickness(gdir): the glacier directory to process """ - base_url = default_base_url - + # Find out which file(s) we need gdf = _get_lookup_thickness() - cp = shpg.Point(gdir.cenlon, gdir.cenlat) sel = gdf.loc[gdf.contains(cp)] if len(sel) > 1: @@ -52,14 +51,27 @@ def add_millan_thickness(gdir): raise InvalidWorkflowError(f'There seems to be no Millan file for this ' f'glacier: {gdir.rgi_id}') - url = base_url + sel['path'].iloc[0] + # Fetch it + url = default_base_url + sel['thickness'].iloc[0] input_file = utils.file_downloader(url) + # Subset to avoid mega files dsb = salem.GeoTiff(input_file) x0, x1, y0, y1 = gdir.grid.extent dsb.set_subset(corners=((x0, y0), (x1, y1)), crs=gdir.grid.proj, margin=5) + # Read the data and prevent bad surprises thick = dsb.get_vardata().astype(np.float64) + # Nans with 0 + thick[~ np.isfinite(thick)] = 0 + nmax = np.nanmax(thick) + if nmax == np.inf: + # Replace inf with 0 + thick[thick == nmax] = 0 + # Replace negative values with 0 + thick[thick < 0] = 0 + + # Reproject now thick = gdir.grid.map_gridded_data(thick, dsb.grid, interp='linear').filled(np.nan) # We mask zero ice as nodata @@ -76,5 +88,4 @@ def add_millan_thickness(gdir): v.units = 'm' ln = 'Ice thickness from Millan et al. 2022' v.long_name = ln - v.base_url = base_url v[:] = thick From ed035ad55aedebf86daef324d187a1ec8bc68046 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Wed, 20 Jul 2022 15:39:59 +0200 Subject: [PATCH 03/11] Also allow for multiple files --- oggm/shop/millan22.py | 58 ++++++++++++++++++++++++------------------- 1 file changed, 32 insertions(+), 26 deletions(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index d23ceb2e2..e52ed9734 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -45,37 +45,43 @@ def add_millan_thickness(gdir): gdf = _get_lookup_thickness() cp = shpg.Point(gdir.cenlon, gdir.cenlat) sel = gdf.loc[gdf.contains(cp)] - if len(sel) > 1: - raise NotImplementedError('Multifile Millan not implemented yet') if len(sel) == 0: raise InvalidWorkflowError(f'There seems to be no Millan file for this ' f'glacier: {gdir.rgi_id}') - # Fetch it - url = default_base_url + sel['thickness'].iloc[0] - input_file = utils.file_downloader(url) - - # Subset to avoid mega files - dsb = salem.GeoTiff(input_file) - x0, x1, y0, y1 = gdir.grid.extent - dsb.set_subset(corners=((x0, y0), (x1, y1)), crs=gdir.grid.proj, margin=5) - - # Read the data and prevent bad surprises - thick = dsb.get_vardata().astype(np.float64) - # Nans with 0 - thick[~ np.isfinite(thick)] = 0 - nmax = np.nanmax(thick) - if nmax == np.inf: - # Replace inf with 0 - thick[thick == nmax] = 0 - # Replace negative values with 0 - thick[thick < 0] = 0 - - # Reproject now - thick = gdir.grid.map_gridded_data(thick, dsb.grid, interp='linear').filled(np.nan) + # We may have more than one file + total_thick = 0 + for i, s in sel.iterrows(): + # Fetch it + url = default_base_url + s['thickness'] + input_file = utils.file_downloader(url) + + # Subset to avoid mega files + dsb = salem.GeoTiff(input_file) + x0, x1, y0, y1 = gdir.grid.extent + dsb.set_subset(corners=((x0, y0), (x1, y1)), crs=gdir.grid.proj, margin=5) + + # Read the data and prevent bad surprises + thick = dsb.get_vardata().astype(np.float64) + # Nans with 0 + thick[~ np.isfinite(thick)] = 0 + nmax = np.nanmax(thick) + if nmax == np.inf: + # Replace inf with 0 + thick[thick == nmax] = 0 + # Replace negative values with 0 + thick[thick < 0] = 0 + + if np.nansum(thick) == 0: + # No need to continue + continue + + # Reproject now + r_thick = gdir.grid.map_gridded_data(thick, dsb.grid, interp='linear') + total_thick += r_thick.filled(0) # We mask zero ice as nodata - thick = np.where(thick == 0, np.NaN, thick) + total_thick = np.where(total_thick == 0, np.NaN, total_thick) # Write with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc: @@ -88,4 +94,4 @@ def add_millan_thickness(gdir): v.units = 'm' ln = 'Ice thickness from Millan et al. 2022' v.long_name = ln - v[:] = thick + v[:] = total_thick From 6b70b906e2a5c85e6d91dc498b3e01cdef8156da Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 21 Jul 2022 16:39:33 +0200 Subject: [PATCH 04/11] Add statistics --- oggm/shop/millan22.py | 66 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index e52ed9734..417876187 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -1,6 +1,9 @@ import logging +import os import numpy as np +import pandas as pd +import xarray as xr import shapely.geometry as shpg try: @@ -13,7 +16,7 @@ except ImportError: pass -from oggm import utils +from oggm import utils, cfg from oggm.exceptions import InvalidWorkflowError # Module logger @@ -95,3 +98,64 @@ def add_millan_thickness(gdir): ln = 'Ice thickness from Millan et al. 2022' v.long_name = ln v[:] = total_thick + +@utils.entity_task(log) +def millan_statistics(gdir): + """Gather statistics about the Millan data interpolated to this glacier. + """ + + d = dict() + + # Easy stats - this should always be possible + d['rgi_id'] = gdir.rgi_id + d['rgi_region'] = gdir.rgi_region + d['rgi_subregion'] = gdir.rgi_subregion + d['rgi_area_km2'] = gdir.rgi_area_km2 + + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: + if 'millan_ice_thickness' not in ds: + d['millan_vol_km3'] = np.NaN + d['millan_area_km2'] = 0 + d['millan_perc_cov'] = 0 + else: + thick = ds['millan_ice_thickness'].where(ds.glacier_mask, np.NaN).load() + d['millan_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9) + d['millan_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) + d['millan_perc_cov'] = float(d['millan_area_km2'] / gdir.rgi_area_km2) + + return d + + +@utils.global_task(log) +def compile_millan_statistics(gdirs, filesuffix='', path=True): + """Gather as much statistics as possible about a list of glaciers. + + It can be used to do result diagnostics and other stuffs. If the data + necessary for a statistic is not available (e.g.: flowlines length) it + will simply be ignored. + + Parameters + ---------- + gdirs : list of :py:class:`oggm.GlacierDirectory` objects + the glacier directories to process + filesuffix : str + add suffix to output file + path : str, bool + Set to "True" in order to store the info in the working directory + Set to a path to store the file to your chosen location + """ + from oggm.workflow import execute_entity_task + + out_df = execute_entity_task(millan_statistics, gdirs) + + out = pd.DataFrame(out_df).set_index('rgi_id') + + if path: + if path is True: + out.to_csv(os.path.join(cfg.PATHS['working_dir'], + ('millan_statistics' + + filesuffix + '.csv'))) + else: + out.to_csv(path) + + return out From 0811b5f1875f2bba7aca574bd16ce0c6c2d584bf Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 21 Jul 2022 17:03:52 +0200 Subject: [PATCH 05/11] Add statistics --- oggm/shop/millan22.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index 417876187..3ccf39d7d 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -112,16 +112,21 @@ def millan_statistics(gdir): d['rgi_subregion'] = gdir.rgi_subregion d['rgi_area_km2'] = gdir.rgi_area_km2 - with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: - if 'millan_ice_thickness' not in ds: - d['millan_vol_km3'] = np.NaN - d['millan_area_km2'] = 0 - d['millan_perc_cov'] = 0 - else: - thick = ds['millan_ice_thickness'].where(ds.glacier_mask, np.NaN).load() - d['millan_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9) - d['millan_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) - d['millan_perc_cov'] = float(d['millan_area_km2'] / gdir.rgi_area_km2) + try: + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: + if 'millan_ice_thickness' not in ds: + d['millan_vol_km3'] = 0 + d['millan_area_km2'] = 0 + d['millan_perc_cov'] = 0 + else: + thick = ds['millan_ice_thickness'].where(ds.glacier_mask, np.NaN).load() + d['millan_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9) + d['millan_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) + d['millan_perc_cov'] = float(d['millan_area_km2'] / gdir.rgi_area_km2) + except FileNotFoundError: + d['millan_vol_km3'] = 0 + d['millan_area_km2'] = 0 + d['millan_perc_cov'] = 0 return d From 78a0bded161b2977aa77e4e1cb4044d397264648 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Fri, 22 Jul 2022 11:17:16 +0200 Subject: [PATCH 06/11] Be more inclusive with errors --- oggm/shop/millan22.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index 3ccf39d7d..806df3e2b 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -111,22 +111,18 @@ def millan_statistics(gdir): d['rgi_region'] = gdir.rgi_region d['rgi_subregion'] = gdir.rgi_subregion d['rgi_area_km2'] = gdir.rgi_area_km2 + d['millan_vol_km3'] = 0 + d['millan_area_km2'] = 0 + d['millan_perc_cov'] = 0 try: with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: - if 'millan_ice_thickness' not in ds: - d['millan_vol_km3'] = 0 - d['millan_area_km2'] = 0 - d['millan_perc_cov'] = 0 - else: - thick = ds['millan_ice_thickness'].where(ds.glacier_mask, np.NaN).load() - d['millan_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9) - d['millan_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) - d['millan_perc_cov'] = float(d['millan_area_km2'] / gdir.rgi_area_km2) - except FileNotFoundError: - d['millan_vol_km3'] = 0 - d['millan_area_km2'] = 0 - d['millan_perc_cov'] = 0 + thick = ds['millan_ice_thickness'].where(ds['glacier_mask'], np.NaN).load() + d['millan_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9) + d['millan_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) + d['millan_perc_cov'] = float(d['millan_area_km2'] / gdir.rgi_area_km2) + except (FileNotFoundError, AttributeError, KeyError): + pass return d From 521bc7c089a78893904324a8fe3c573ec9cd8a88 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Fri, 2 Sep 2022 15:14:00 +0200 Subject: [PATCH 07/11] Also velocity --- oggm/shop/millan22.py | 230 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 199 insertions(+), 31 deletions(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index 806df3e2b..2fd24f460 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -25,38 +25,54 @@ default_base_url = 'https://cluster.klima.uni-bremen.de/~oggm/velocities/millan22/' _lookup_thickness = None +_lookup_velocity = None + def _get_lookup_thickness(): global _lookup_thickness if _lookup_thickness is None: - fname = default_base_url + 'millan22_thickness_lookup_shp.zip' + fname = default_base_url + 'millan22_thickness_lookup_shp_20220902.zip' _lookup_thickness = gpd.read_file('zip://' + utils.file_downloader(fname)) return _lookup_thickness -@utils.entity_task(log, writes=['gridded_data']) -def add_millan_thickness(gdir): - """Add the Millan 22 thickness data to this glacier directory. +def _get_lookup_velocity(): + global _lookup_velocity + if _lookup_velocity is None: + fname = default_base_url + 'millan22_velocity_lookup_shp_20220902.zip' + _lookup_velocity = gpd.read_file('zip://' + utils.file_downloader(fname)) + return _lookup_velocity + + +def _filter(ds): + # Read the data and prevent bad surprises + data = ds.get_vardata().astype(np.float64) + # Nans with 0 + data[~ np.isfinite(data)] = 0 + nmax = np.nanmax(data) + if nmax == np.inf: + # Replace inf with 0 + data[data == nmax] = 0 + + return data + + +def _filter_and_reproj(gdir, var, gdf, allow_neg=True): + """ Same code for thickness and error Parameters ---------- - gdir ::py:class:`oggm.GlacierDirectory` - the glacier directory to process + var : 'thickness' or 'err' + gdf : the lookup """ - # Find out which file(s) we need - gdf = _get_lookup_thickness() - cp = shpg.Point(gdir.cenlon, gdir.cenlat) - sel = gdf.loc[gdf.contains(cp)] - if len(sel) == 0: - raise InvalidWorkflowError(f'There seems to be no Millan file for this ' - f'glacier: {gdir.rgi_id}') - # We may have more than one file - total_thick = 0 - for i, s in sel.iterrows(): + total_data = 0 + grids_used = [] + files_used = [] + for i, s in gdf.iterrows(): # Fetch it - url = default_base_url + s['thickness'] + url = default_base_url + s[var] input_file = utils.file_downloader(url) # Subset to avoid mega files @@ -64,28 +80,56 @@ def add_millan_thickness(gdir): x0, x1, y0, y1 = gdir.grid.extent dsb.set_subset(corners=((x0, y0), (x1, y1)), crs=gdir.grid.proj, margin=5) - # Read the data and prevent bad surprises - thick = dsb.get_vardata().astype(np.float64) - # Nans with 0 - thick[~ np.isfinite(thick)] = 0 - nmax = np.nanmax(thick) - if nmax == np.inf: - # Replace inf with 0 - thick[thick == nmax] = 0 - # Replace negative values with 0 - thick[thick < 0] = 0 - - if np.nansum(thick) == 0: + data = _filter(dsb) + if not allow_neg: + # Replace negative values with 0 + data[data < 0] = 0 + + if np.nansum(data) == 0: # No need to continue continue # Reproject now - r_thick = gdir.grid.map_gridded_data(thick, dsb.grid, interp='linear') - total_thick += r_thick.filled(0) + r_data = gdir.grid.map_gridded_data(data, dsb.grid, interp='linear') + total_data += r_data.filled(0) + grids_used.append(dsb) + files_used.append(s.file_id) + + return total_data, files_used, grids_used + + +@utils.entity_task(log, writes=['gridded_data']) +def add_millan_thickness(gdir, add_error=False): + """Add the Millan 22 thickness data to this glacier directory. + + Parameters + ---------- + gdir : :py:class:`oggm.GlacierDirectory` + the glacier directory to process + add_error : bool + add the error data or not + """ + + # Find out which file(s) we need + gdf = _get_lookup_thickness() + cp = shpg.Point(gdir.cenlon, gdir.cenlat) + sel = gdf.loc[gdf.contains(cp)] + if len(sel) == 0: + raise InvalidWorkflowError(f'There seems to be no Millan file for this ' + f'glacier: {gdir.rgi_id}') + + total_thick, files_used, _ = _filter_and_reproj(gdir, 'thickness', sel, + allow_neg=False) # We mask zero ice as nodata total_thick = np.where(total_thick == 0, np.NaN, total_thick) + if add_error: + total_err, _, _ = _filter_and_reproj(gdir, 'err', sel, allow_neg=False) + total_err[~ np.isfinite(total_thick)] = np.NaN + # Error cannot be larger than ice thickness itself + total_err = utils.clip_max(total_err, total_thick) + # Write with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc: @@ -97,8 +141,119 @@ def add_millan_thickness(gdir): v.units = 'm' ln = 'Ice thickness from Millan et al. 2022' v.long_name = ln + data_str = ' '.join(files_used) if len(files_used) > 1 else files_used[0] + v.data_source = data_str v[:] = total_thick + if add_error: + vn = 'millan_ice_thickness_err' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x',), zlib=True) + v.units = 'm' + ln = 'Ice thickness error from Millan et al. 2022' + v.long_name = ln + v.data_source = data_str + v[:] = total_err + + +@utils.entity_task(log, writes=['gridded_data']) +def add_millan_velocity(gdir, add_error=False): + """Add the Millan 22 thickness data to this glacier directory. + + Parameters + ---------- + gdir : :py:class:`oggm.GlacierDirectory` + the glacier directory to process + add_error : bool + add the error data or not + """ + + # Find out which file(s) we need + gdf = _get_lookup_velocity() + cp = shpg.Point(gdir.cenlon, gdir.cenlat) + sel = gdf.loc[gdf.contains(cp)] + if len(sel) == 0: + raise InvalidWorkflowError(f'There seems to be no Millan file for this ' + f'glacier: {gdir.rgi_id}') + + vel, files, grids = _filter_and_reproj(gdir, 'v', sel, allow_neg=False) + assert len(grids) == 1, 'Multiple velocity grids - dont know what to do.' + sel = sel.loc[sel.file_id == files[0]] + vx, _, gridsx = _filter_and_reproj(gdir, 'vx', sel) + vy, _, gridsy = _filter_and_reproj(gdir, 'vy', sel) + + dsx = gridsx[0] + dsy = gridsy[0] + grid_vel = dsx.grid + proj_vel = grid_vel.proj + grid_gla = gdir.grid + + # Get the coords at t0 + xx0, yy0 = grid_vel.center_grid.xy_coordinates + + # Compute coords at t1 + xx1 = _filter(dsx) + yy1 = _filter(dsy) + xx1 += xx0 + yy1 += yy0 + + # Transform both to glacier proj + xx0, yy0 = salem.transform_proj(proj_vel, grid_gla.proj, xx0, yy0) + xx1, yy1 = salem.transform_proj(proj_vel, grid_gla.proj, xx1, yy1) + + # Compute velocities from there + vx = xx1 - xx0 + vy = yy1 - yy0 + + # And transform to local map + vx = grid_gla.map_gridded_data(vx, grid=grid_vel, interp='linear') + vy = grid_gla.map_gridded_data(vy, grid=grid_vel, interp='linear') + + # Scale back to match velocity + new_vel = np.sqrt(vx**2 + vy**2) + p_ok = np.isfinite(new_vel) & (new_vel > 1e-5) # avoid div by zero + vx[p_ok] = vx[p_ok] * vel[p_ok] / new_vel[p_ok] + vy[p_ok] = vy[p_ok] * vel[p_ok] / new_vel[p_ok] + + # Write + with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc: + + vn = 'millan_v' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x', ), zlib=True) + v.units = 'm' + ln = 'Ice velocity from Millan et al. 2022' + v.long_name = ln + v.data_source = files[0] + v[:] = vel + + vn = 'millan_vx' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x', ), zlib=True) + v.units = 'm' + ln = 'Ice velocity in map x direction from Millan et al. 2022' + v.long_name = ln + v.data_source = files[0] + v[:] = vx + + vn = 'millan_vy' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x', ), zlib=True) + v.units = 'm' + ln = 'Ice velocity in map y direction from Millan et al. 2022' + v.long_name = ln + v.data_source = files[0] + v[:] = vy + + @utils.entity_task(log) def millan_statistics(gdir): """Gather statistics about the Millan data interpolated to this glacier. @@ -121,6 +276,19 @@ def millan_statistics(gdir): d['millan_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9) d['millan_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) d['millan_perc_cov'] = float(d['millan_area_km2'] / gdir.rgi_area_km2) + + if 'millan_ice_thickness_err' in ds: + err = ds['millan_ice_thickness_err'].where(ds['glacier_mask'], np.NaN).load() + d['millan_vol_err_km3'] = float(err.sum() * gdir.grid.dx ** 2 * 1e-9) + except (FileNotFoundError, AttributeError, KeyError): + pass + + try: + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: + v = ds['millan_v'].where(ds['glacier_mask'], np.NaN).load() + d['millan_avg_vel'] = np.nanmean(v) + d['millan_max_vel'] = np.nanmax(v) + except (FileNotFoundError, AttributeError, KeyError): pass From 3ed364d5057840c34a5a6b2b6b5bec0ff84d8019 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Fri, 2 Sep 2022 17:55:24 +0200 Subject: [PATCH 08/11] Add v error + fixes --- oggm/shop/its_live.py | 3 +-- oggm/shop/millan22.py | 55 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 52 insertions(+), 6 deletions(-) diff --git a/oggm/shop/its_live.py b/oggm/shop/its_live.py index 63e1e2998..1c71abfec 100644 --- a/oggm/shop/its_live.py +++ b/oggm/shop/its_live.py @@ -78,7 +78,6 @@ def find_region(gdir): def _reproject_and_scale(gdir, do_error=False): """Reproject and scale itslive data, avoid code duplication for error""" - reg = find_region(gdir) if reg is None: raise InvalidWorkflowError('There does not seem to be its_live data ' @@ -142,7 +141,7 @@ def _reproject_and_scale(gdir, do_error=False): # Scale back velocities - https://github.com/OGGM/oggm/issues/1014 new_vel = np.sqrt(vx**2 + vy**2) - p_ok = new_vel > 1e-5 # avoid div by zero + p_ok = new_vel > 1 # avoid div by zero vx[p_ok] = vx[p_ok] * orig_vel[p_ok] / new_vel[p_ok] vy[p_ok] = vy[p_ok] * orig_vel[p_ok] / new_vel[p_ok] diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index 2fd24f460..1cbb0a71d 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -1,5 +1,6 @@ import logging import os +import warnings import numpy as np import pandas as pd @@ -170,6 +171,11 @@ def add_millan_velocity(gdir, add_error=False): add the error data or not """ + if gdir.rgi_region in ['05']: + raise NotImplementedError('Millan 22 does not provide velocity ' + 'products for Greenland - we would need to ' + 'implement MEASURES in the shop for this.') + # Find out which file(s) we need gdf = _get_lookup_velocity() cp = shpg.Point(gdir.cenlon, gdir.cenlat) @@ -212,10 +218,22 @@ def add_millan_velocity(gdir, add_error=False): vy = grid_gla.map_gridded_data(vy, grid=grid_vel, interp='linear') # Scale back to match velocity - new_vel = np.sqrt(vx**2 + vy**2) - p_ok = np.isfinite(new_vel) & (new_vel > 1e-5) # avoid div by zero - vx[p_ok] = vx[p_ok] * vel[p_ok] / new_vel[p_ok] - vy[p_ok] = vy[p_ok] * vel[p_ok] / new_vel[p_ok] + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning) + new_vel = np.sqrt(vx**2 + vy**2) + p_ok = np.isfinite(new_vel) & (new_vel > 1) # avoid div by zero + scaler = vel[p_ok] / new_vel[p_ok] + vx[p_ok] = vx[p_ok] * scaler + vy[p_ok] = vy[p_ok] * scaler + + vx = vx.filled(np.NaN) + vy = vy.filled(np.NaN) + + if add_error: + err_vx, _, _ = _filter_and_reproj(gdir, 'err_vx', sel, allow_neg=False) + err_vy, _, _ = _filter_and_reproj(gdir, 'err_vy', sel, allow_neg=False) + err_vx[p_ok] = err_vx[p_ok] * scaler + err_vy[p_ok] = err_vy[p_ok] * scaler # Write with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc: @@ -253,6 +271,29 @@ def add_millan_velocity(gdir, add_error=False): v.data_source = files[0] v[:] = vy + if add_error: + vn = 'millan_err_vx' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x',), zlib=True) + v.units = 'm' + ln = 'Ice velocity error in map x direction from Millan et al. 2022' + v.long_name = ln + v.data_source = files[0] + v[:] = err_vx + + vn = 'millan_err_vy' + if vn in nc.variables: + v = nc.variables[vn] + else: + v = nc.createVariable(vn, 'f4', ('y', 'x',), zlib=True) + v.units = 'm' + ln = 'Ice velocity error in map y direction from Millan et al. 2022' + v.long_name = ln + v.data_source = files[0] + v[:] = err_vy + @utils.entity_task(log) def millan_statistics(gdir): @@ -289,6 +330,12 @@ def millan_statistics(gdir): d['millan_avg_vel'] = np.nanmean(v) d['millan_max_vel'] = np.nanmax(v) + if 'millan_err_vx' in ds: + err_vx = ds['millan_err_vx'].where(ds['glacier_mask'], np.NaN).load() + err_vy = ds['millan_err_vy'].where(ds['glacier_mask'], np.NaN).load() + err = (err_vx**2 + err_vy**2)**0.5 + d['millan_avg_err_vel'] = np.nanmean(err) + except (FileNotFoundError, AttributeError, KeyError): pass From e4ae01121e5d5bbfc011b0e67ad11590a2cd9002 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sun, 4 Sep 2022 11:32:24 +0200 Subject: [PATCH 09/11] Add tests --- oggm/shop/millan22.py | 4 +- oggm/tests/conftest.py | 1 + oggm/tests/test_shop.py | 120 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 122 insertions(+), 3 deletions(-) diff --git a/oggm/shop/millan22.py b/oggm/shop/millan22.py index 1cbb0a71d..5622d4c50 100644 --- a/oggm/shop/millan22.py +++ b/oggm/shop/millan22.py @@ -100,7 +100,7 @@ def _filter_and_reproj(gdir, var, gdf, allow_neg=True): @utils.entity_task(log, writes=['gridded_data']) -def add_millan_thickness(gdir, add_error=False): +def thickness_to_gdir(gdir, add_error=False): """Add the Millan 22 thickness data to this glacier directory. Parameters @@ -160,7 +160,7 @@ def add_millan_thickness(gdir, add_error=False): @utils.entity_task(log, writes=['gridded_data']) -def add_millan_velocity(gdir, add_error=False): +def velocity_to_gdir(gdir, add_error=False): """Add the Millan 22 thickness data to this glacier directory. Parameters diff --git a/oggm/tests/conftest.py b/oggm/tests/conftest.py index 5d4b568f6..6f4d2288d 100644 --- a/oggm/tests/conftest.py +++ b/oggm/tests/conftest.py @@ -155,6 +155,7 @@ def secure_url_retrieve(url, *args, **kwargs): 'cluster.klima.uni-bremen.de/~oggm/test_gdirs/' in url or 'cluster.klima.uni-bremen.de/~oggm/demo_gdirs/' in url or 'cluster.klima.uni-bremen.de/~oggm/test_climate/' in url or + 'cluster.klima.uni-bremen.de/~oggm/test_files/' in url or 'klima.uni-bremen.de/~oggm/climate/cru/cru_cl2.nc.zip' in url or 'klima.uni-bremen.de/~oggm/geodetic_ref_mb' in url or base_extra_L4 in url or diff --git a/oggm/tests/test_shop.py b/oggm/tests/test_shop.py index a8ea39475..ff6838828 100644 --- a/oggm/tests/test_shop.py +++ b/oggm/tests/test_shop.py @@ -12,7 +12,7 @@ import pandas as pd from oggm import utils from oggm.utils import get_demo_file -from oggm.shop import its_live, rgitopo, bedtopo +from oggm.shop import its_live, rgitopo, bedtopo, millan22 from oggm.core import gis, centerlines, massbalance from oggm import cfg, tasks, workflow @@ -107,6 +107,124 @@ def test_repro_to_glacier(self, class_case_dir, monkeypatch): plt.show() +class Test_millan22: + + @pytest.mark.slow + def test_thick_to_glacier(self, class_case_dir, monkeypatch): + + # Init + cfg.initialize() + cfg.PATHS['working_dir'] = class_case_dir + cfg.PARAMS['use_intersects'] = False + cfg.PATHS['dem_file'] = get_demo_file('dem_Columbia.tif') + cfg.PARAMS['border'] = 10 + + entity = gpd.read_file(get_demo_file('RGI60-01.10689.shp')).iloc[0] + gdir = oggm.GlacierDirectory(entity) + tasks.define_glacier_region(gdir) + tasks.glacier_masks(gdir) + + # use our files + base_url = 'https://cluster.klima.uni-bremen.de/~oggm/test_files/millan22/' + monkeypatch.setattr(millan22, 'default_base_url', base_url) + + millan22.thickness_to_gdir(gdir) + + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: + mask = ds.glacier_mask.data.astype(bool) + thick = ds.millan_ice_thickness.where(mask).data + + # Simply some coverage and sanity checks + assert np.isfinite(thick).sum() / mask.sum() > 0.98 + assert np.nanmax(thick) > 800 + assert np.nansum(thick) * gdir.grid.dx**2 * 1e-9 > 174 + + df = millan22.compile_millan_statistics([gdir]).iloc[0] + assert df['millan_perc_cov'] > 0.95 + assert df['millan_vol_km3'] > 174 + + if DO_PLOT: + import matplotlib.pyplot as plt + + smap = salem.Map(gdir.grid.center_grid, countries=False) + smap.set_shapefile(gdir.read_shapefile('outlines')) + + with warnings.catch_warnings(): + warnings.filterwarnings('ignore', category=RuntimeWarning) + smap.set_topography(gdir.get_filepath('dem')) + + smap.set_data(thick) + smap.set_plot_params(cmap='Reds', vmin=None, vmax=None) + + f, ax = plt.subplots() + smap.visualize(ax=ax, title='M22 thickness', + cbar_title='m') + plt.show() + + @pytest.mark.slow + def test_vel_to_glacier(self, class_case_dir, monkeypatch): + + # Init + cfg.initialize() + cfg.PATHS['working_dir'] = class_case_dir + cfg.PARAMS['use_intersects'] = False + cfg.PATHS['dem_file'] = get_demo_file('dem_Columbia.tif') + cfg.PARAMS['border'] = 10 + + entity = gpd.read_file(get_demo_file('RGI60-01.10689.shp')).iloc[0] + gdir = oggm.GlacierDirectory(entity) + tasks.define_glacier_region(gdir) + tasks.glacier_masks(gdir) + + # use our files + base_url = 'https://cluster.klima.uni-bremen.de/~oggm/test_files/millan22/' + monkeypatch.setattr(millan22, 'default_base_url', base_url) + + millan22.velocity_to_gdir(gdir) + + with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: + mask = ds.glacier_mask.data.astype(bool) + v = ds.millan_v.where(mask).data + vx = ds.millan_vx.where(mask).data + vy = ds.millan_vy.where(mask).data + + # Simply some coverage and sanity checks + assert np.isfinite(v).sum() / mask.sum() > 0.98 + assert np.nanmax(v) > 2000 + + df = millan22.compile_millan_statistics([gdir]).iloc[0] + assert df['millan_avg_vel'] > 180 + assert df['millan_max_vel'] > 2000 + + if DO_PLOT: + import matplotlib.pyplot as plt + + smap = salem.Map(gdir.grid.center_grid, countries=False) + smap.set_shapefile(gdir.read_shapefile('outlines')) + + with warnings.catch_warnings(): + warnings.filterwarnings('ignore', category=RuntimeWarning) + smap.set_topography(gdir.get_filepath('dem')) + + smap.set_data(v) + smap.set_plot_params(cmap='Blues', vmin=None, vmax=None) + + xx, yy = gdir.grid.center_grid.xy_coordinates + xx, yy = smap.grid.transform(xx, yy, crs=gdir.grid.proj) + + yy = yy[2::5, 2::5] + xx = xx[2::5, 2::5] + vx = vx[2::5, 2::5] + vy = vy[2::5, 2::5] + + f, ax = plt.subplots() + smap.visualize(ax=ax, title='M22 velocity', + cbar_title='m yr-1') + ax.quiver(xx, yy, vx, vy) + plt.show() + plt.show() + + class Test_rgitopo: def test_from_dem(self, class_case_dir, monkeypatch): From 2858145e88f8fa21093674856d0675a60014cf04 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sun, 4 Sep 2022 11:33:49 +0200 Subject: [PATCH 10/11] whatsnew --- docs/whats-new.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/whats-new.rst b/docs/whats-new.rst index 28f96ae8c..fdb4febcb 100644 --- a/docs/whats-new.rst +++ b/docs/whats-new.rst @@ -40,6 +40,8 @@ Enhancements Level 5 now replaces level 4 and creates the minigdirs (where only the files for a model run are kept and no inversion is possible anymore) (:pull:`1425`). By `Patrick Schmitt `_ +- Added support for Milland 22 velocity and thickness in the shop (:pull:`1443`). + By `Fabien Maussion `_ Bug fixes ~~~~~~~~~ From 2357810e06f639904a592b9a0812c2ba823ff51e Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Mon, 5 Sep 2022 09:36:59 +0200 Subject: [PATCH 11/11] update test --- oggm/tests/test_shop.py | 75 +++-------------------------------------- 1 file changed, 5 insertions(+), 70 deletions(-) diff --git a/oggm/tests/test_shop.py b/oggm/tests/test_shop.py index ff6838828..eb9747eff 100644 --- a/oggm/tests/test_shop.py +++ b/oggm/tests/test_shop.py @@ -110,7 +110,7 @@ def test_repro_to_glacier(self, class_case_dir, monkeypatch): class Test_millan22: @pytest.mark.slow - def test_thick_to_glacier(self, class_case_dir, monkeypatch): + def test_thickvel_to_glacier(self, class_case_dir, monkeypatch): # Init cfg.initialize() @@ -139,47 +139,7 @@ def test_thick_to_glacier(self, class_case_dir, monkeypatch): assert np.nanmax(thick) > 800 assert np.nansum(thick) * gdir.grid.dx**2 * 1e-9 > 174 - df = millan22.compile_millan_statistics([gdir]).iloc[0] - assert df['millan_perc_cov'] > 0.95 - assert df['millan_vol_km3'] > 174 - - if DO_PLOT: - import matplotlib.pyplot as plt - - smap = salem.Map(gdir.grid.center_grid, countries=False) - smap.set_shapefile(gdir.read_shapefile('outlines')) - - with warnings.catch_warnings(): - warnings.filterwarnings('ignore', category=RuntimeWarning) - smap.set_topography(gdir.get_filepath('dem')) - - smap.set_data(thick) - smap.set_plot_params(cmap='Reds', vmin=None, vmax=None) - - f, ax = plt.subplots() - smap.visualize(ax=ax, title='M22 thickness', - cbar_title='m') - plt.show() - - @pytest.mark.slow - def test_vel_to_glacier(self, class_case_dir, monkeypatch): - - # Init - cfg.initialize() - cfg.PATHS['working_dir'] = class_case_dir - cfg.PARAMS['use_intersects'] = False - cfg.PATHS['dem_file'] = get_demo_file('dem_Columbia.tif') - cfg.PARAMS['border'] = 10 - - entity = gpd.read_file(get_demo_file('RGI60-01.10689.shp')).iloc[0] - gdir = oggm.GlacierDirectory(entity) - tasks.define_glacier_region(gdir) - tasks.glacier_masks(gdir) - - # use our files - base_url = 'https://cluster.klima.uni-bremen.de/~oggm/test_files/millan22/' - monkeypatch.setattr(millan22, 'default_base_url', base_url) - + # Velocity millan22.velocity_to_gdir(gdir) with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds: @@ -192,37 +152,12 @@ def test_vel_to_glacier(self, class_case_dir, monkeypatch): assert np.isfinite(v).sum() / mask.sum() > 0.98 assert np.nanmax(v) > 2000 + # Stats df = millan22.compile_millan_statistics([gdir]).iloc[0] assert df['millan_avg_vel'] > 180 assert df['millan_max_vel'] > 2000 - - if DO_PLOT: - import matplotlib.pyplot as plt - - smap = salem.Map(gdir.grid.center_grid, countries=False) - smap.set_shapefile(gdir.read_shapefile('outlines')) - - with warnings.catch_warnings(): - warnings.filterwarnings('ignore', category=RuntimeWarning) - smap.set_topography(gdir.get_filepath('dem')) - - smap.set_data(v) - smap.set_plot_params(cmap='Blues', vmin=None, vmax=None) - - xx, yy = gdir.grid.center_grid.xy_coordinates - xx, yy = smap.grid.transform(xx, yy, crs=gdir.grid.proj) - - yy = yy[2::5, 2::5] - xx = xx[2::5, 2::5] - vx = vx[2::5, 2::5] - vy = vy[2::5, 2::5] - - f, ax = plt.subplots() - smap.visualize(ax=ax, title='M22 velocity', - cbar_title='m yr-1') - ax.quiver(xx, yy, vx, vy) - plt.show() - plt.show() + assert df['millan_perc_cov'] > 0.95 + assert df['millan_vol_km3'] > 174 class Test_rgitopo: