Add hugonnet dhdt to shop and both hugonnet and millan to prepro leve…

…ls (#1529)

docs/whats-new.rst

@@ -55,6 +55,8 @@ Enhancements
   By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 - Added support for Millan et al 2022 velocity and thickness in the shop (:pull:`1443`).
   By `Fabien Maussion <https://github.com/fmaussion>`_
+- Added support for Hugonnet et al 2021 dhdt in the shop (:pull:`1529`).
+  By `Fabien Maussion <https://github.com/fmaussion>`_
 - Added trapezoidal downstream line (:pull:`1491`). Can be selected with
   ``cfg.PARAMS['downstream_line_shape']``, with the options ``'parabol'`` (default)
   or ``'trapezoidal'`` before calling ``init_present_time_glacier(gdir)``.

oggm/cli/prepro_levels.py

@@ -81,8 +81,10 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
                       match_geodetic_mb_per_glacier=False,
                       evolution_model='fl_sia',
                       centerlines_only=False, override_params=None,
-                      add_consensus=False, start_level=None,
-                      start_base_url=None, max_level=5, ref_tstars_base_url='',
+                      add_consensus_thickness=False, add_millan_thickness=False,
+                      add_millan_velocity=False, add_hugonnet_dhdt=False,
+                      start_level=None, start_base_url=None, max_level=5,
+                      ref_tstars_base_url='',
                       logging_level='WORKFLOW', disable_dl_verify=False,
                       dynamic_spinup=False, err_dmdtda_scaling_factor=1,
                       dynamic_spinup_start_year=1979,
@@ -139,9 +141,18 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
     evolution_model : str
         which geometry evolution model to use: `fl_sia` (default),
         or `massredis` (mass redistribution curve).
-    add_consensus : bool
+    add_consensus_thickness : bool
         adds (reprojects) the consensus estimates thickness to the glacier
         directories. With elev_bands=True, the data will also be binned.
+    add_millan_thickness : bool
+        adds (reprojects) the millan thickness to the glacier
+        directories. With elev_bands=True, the data will also be binned.
+    add_millan_velocity : bool
+        adds (reprojects) the millan velocity to the glacier
+        directories. With elev_bands=True, the data will also be binned.
+    add_hugonnet_dhdt : bool
+        adds (reprojects) the hugonnet dhdt maps to the glacier
+        directories. With elev_bands=True, the data will also be binned.
     start_level : int
         the pre-processed level to start from (default is to start from
         scratch). If set, you'll need to indicate start_base_url as well.
@@ -430,17 +441,33 @@ def _time_log():
         # Centerlines OGGM
         workflow.execute_entity_task(tasks.glacier_masks, gdirs_cent)
 
-        if add_consensus:
+        bin_variables = []
+        if add_consensus_thickness:
             from oggm.shop.bedtopo import add_consensus_thickness
             workflow.execute_entity_task(add_consensus_thickness, gdirs_band)
             workflow.execute_entity_task(add_consensus_thickness, gdirs_cent)
-
-            # Elev bands with var data
-            vn = 'consensus_ice_thickness'
+            bin_variables.append('consensus_ice_thickness')
+        if add_millan_thickness:
+            from oggm.shop.millan22 import thickness_to_gdir
+            workflow.execute_entity_task(thickness_to_gdir, gdirs_band)
+            workflow.execute_entity_task(thickness_to_gdir, gdirs_cent)
+            bin_variables.append('millan_ice_thickness')
+        if add_millan_velocity:
+            from oggm.shop.millan22 import velocity_to_gdir
+            workflow.execute_entity_task(velocity_to_gdir, gdirs_band)
+            workflow.execute_entity_task(velocity_to_gdir, gdirs_cent)
+            bin_variables.append('millan_v')
+        if add_hugonnet_dhdt:
+            from oggm.shop.hugonnet_maps import hugonnet_to_gdir
+            workflow.execute_entity_task(hugonnet_to_gdir, gdirs_band)
+            workflow.execute_entity_task(hugonnet_to_gdir, gdirs_cent)
+            bin_variables.append('hugonnet_dhdt')
+
+        if bin_variables and gdirs_band:
             workflow.execute_entity_task(tasks.elevation_band_flowline,
-                                         gdirs_band, bin_variables=vn)
+                                         gdirs_band, bin_variables=bin_variables)
             workflow.execute_entity_task(tasks.fixed_dx_elevation_band_flowline,
-                                         gdirs_band, bin_variables=vn)
+                                         gdirs_band, bin_variables=bin_variables)
         else:
             # HH2015 method without it
             task_list = [
@@ -480,6 +507,19 @@ def _time_log():
         opath = os.path.join(sum_dir, 'glacier_statistics_{}.csv'.format(rgi_reg))
         utils.compile_glacier_statistics(gdirs, path=opath)
 
+        if add_millan_thickness or add_millan_velocity:
+            from oggm.shop.millan22 import compile_millan_statistics
+            opath = os.path.join(sum_dir, 'millan_statistics_{}.csv'.format(rgi_reg))
+            compile_millan_statistics(gdirs, path=opath)
+        if add_hugonnet_dhdt:
+            from oggm.shop.hugonnet_maps import compile_hugonnet_statistics
+            opath = os.path.join(sum_dir, 'hugonnet_statistics_{}.csv'.format(rgi_reg))
+            compile_hugonnet_statistics(gdirs, path=opath)
+        if add_consensus_thickness:
+            from oggm.shop.bedtopo import compile_consensus_statistics
+            opath = os.path.join(sum_dir, 'consensus_statistics_{}.csv'.format(rgi_reg))
+            compile_consensus_statistics(gdirs, path=opath)
+
         # And for level 2: shapes
         if len(gdirs_cent) > 0:
             opath = os.path.join(sum_dir, f'centerlines_{rgi_reg}.shp')
@@ -814,9 +854,24 @@ def parse_args(args):
                              'compatible with level 1 folders, after which '
                              'the processing will stop. The default is to use '
                              'the default OGGM DEM.')
-    parser.add_argument('--add-consensus', nargs='?', const=True, default=False,
-                        help='adds (reprojects) the consensus estimates '
-                             'thickness to the glacier directories. '
+    parser.add_argument('--add-consensus-thickness', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the consensus thickness '
+                             'estimates to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-millan-thickness', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the millan thickness '
+                             'estimates to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-millan-velocity', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the millan velocity '
+                             'estimates to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-hugonnet-dhdt', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the millan dhdt '
+                             'maps to the glacier directories. '
                              'With --elev-bands, the data will also be '
                              'binned.')
     parser.add_argument('--demo', nargs='?', const=True, default=False,
@@ -899,7 +954,10 @@ def parse_args(args):
                 centerlines_only=args.centerlines_only,
                 match_regional_geodetic_mb=args.match_regional_geodetic_mb,
                 match_geodetic_mb_per_glacier=args.match_geodetic_mb_per_glacier,
-                add_consensus=args.add_consensus,
+                add_consensus_thickness=args.add_consensus_thickness,
+                add_millan_thickness=args.add_millan_thickness,
+                add_millan_velocity=args.add_millan_velocity,
+                add_hugonnet_dhdt=args.add_hugonnet_dhdt,
                 disable_dl_verify=args.disable_dl_verify,
                 ref_tstars_base_url=args.ref_tstars_base_url,
                 evolution_model=args.evolution_model,

oggm/shop/bedtopo.py

@@ -1,13 +1,16 @@
 import logging
+import os
 
 import numpy as np
+import pandas as pd
+import xarray as xr
 
 try:
     import salem
 except ImportError:
     pass
 
-from oggm import utils
+from oggm import utils, cfg
 
 # Module logger
 log = logging.getLogger(__name__)
@@ -68,3 +71,62 @@ def add_consensus_thickness(gdir, base_url=None):
         v.long_name = ln
         v.base_url = base_url
         v[:] = thick
+
+
+@utils.entity_task(log)
+def consensus_statistics(gdir):
+    """Gather statistics about the consensus 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
+    d['consensus_vol_km3'] = 0
+    d['consensus_area_km2'] = 0
+    d['consensus_perc_cov'] = 0
+
+    try:
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+            thick = ds['consensus_ice_thickness'].where(ds['glacier_mask'], np.NaN).load()
+            d['consensus_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9)
+            d['consensus_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6)
+            d['consensus_perc_cov'] = float(d['consensus_area_km2'] / gdir.rgi_area_km2)
+    except (FileNotFoundError, AttributeError, KeyError):
+        pass
+
+    return d
+
+
+@utils.global_task(log)
+def compile_consensus_statistics(gdirs, filesuffix='', path=True):
+    """Gather as much statistics as possible about a list of glaciers.
+
+    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(consensus_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'],
+                                    ('consensus_statistics' +
+                                     filesuffix + '.csv')))
+        else:
+            out.to_csv(path)
+
+    return out