Refactor distrib funcs and add them to OGGM

.github/workflows/run-tests.yml

@@ -42,12 +42,12 @@ jobs:
       - name: Fix Git-Permissions
         run: git config --global --add safe.directory "$GITHUB_WORKSPACE"
       - name: Test
-        run: ./ci/run_tests.sh ${{ matrix.test-env }}
+        run: ./ci/run_tests.sh "${{ matrix.test-env }}" "${{ matrix.container }}"
       - name: Upload pytest-mpl artifacts
         if: failure()
         uses: actions/upload-artifact@v3
         with:
-          name: pytest-mpl-results-${{ matrix.container }}-${{ matrix.test-env }}
+          name: pytest-mpl-results
           path: /tmp/oggm-mpl-results/
       - name: Upload Coverage
         if: ${{ !contains(matrix.container, 'py3') }}

ci/run_tests.sh

@@ -8,6 +8,8 @@ export OGGM_TEST_ENV="$1"
 export OGGM_TEST_MODE=unset
 export MPLBACKEND=agg
 
+MPL_OUTPUT_OGGM_SUBDIR="$2"
+
 if [[ $OGGM_TEST_ENV == *-* ]]; then
     export OGGM_TEST_MODE="${OGGM_TEST_ENV/*-/}"
     export OGGM_TEST_ENV="${OGGM_TEST_ENV/-*/}"
@@ -34,7 +36,7 @@ export COVERAGE_RCFILE="$PWD/.coveragerc"
 coverage erase
 
 coverage run --source=./oggm --parallel-mode --module \
-    pytest --verbose --mpl-results-path=/tmp/oggm-mpl-results $OGGM_MPL --run-slow --run-test-env $OGGM_TEST_ENV oggm
+    pytest --verbose --mpl-results-path="/tmp/oggm-mpl-results/${MPL_OUTPUT_OGGM_SUBDIR/:/_}/${OGGM_TEST_ENV/:/_}" $OGGM_MPL --run-slow --run-test-env $OGGM_TEST_ENV oggm
 
 coverage combine
 coverage xml

docs/whats-new.rst

@@ -15,6 +15,15 @@ Bug fixes
   special data variables (e.g. ``is_fixed_geometry_spinup``) (:pull:`1563`).
   By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 
+Enhancements
+~~~~~~~~~~~~
+
+- There is now the possibility to compute distributed area and thickness
+  changes from the flowline projections (:pull:`1576`, :pull:`1585`). The
+  functionality is currently in the sandbox and needs to be documented.
+  By `Anouk Vlug <https://github.com/anoukvlug>`_ and
+  `Fabien Maussion <https://github.com/fmaussion>`_
+
 
 v1.6.0 (March 10, 2023)
 -----------------------

oggm/core/flowline.py

@@ -985,12 +985,13 @@ def run_until_and_store(self, y1,
             glacier wide diagnostic files - all other outputs are set
             to constants during "spinup"
         dynamic_spinup_min_ice_thick : float or None
-            if set to an float, additional variables are saved which are useful
+            if set to a float, additional variables are saved which are useful
             in combination with the dynamic spinup. In particular only grid
             points with a minimum ice thickness are considered for the total
             area or the total volume. This is useful to smooth out yearly
             fluctuations when matching to observations. The names of this new
             variables include the suffix _min_h (e.g. 'area_m2_min_h')
+
         Returns
         -------
         geom_ds : xarray.Dataset or None
@@ -1173,6 +1174,8 @@ def run_until_and_store(self, y1,
                 ds.attrs['mb_model_class'] = self.mb_model.__class__.__name__
                 for k, v in self.mb_model.__dict__.items():
                     if np.isscalar(v) and not k.startswith('_'):
+                        if type(v) is bool:
+                            v = str(v)
                         ds.attrs['mb_model_{}'.format(k)] = v
 
                 # Coordinates

oggm/core/inversion.py

@@ -723,20 +723,25 @@ def compute_inversion_velocities(gdir, glen_a=None, fs=None, filesuffix='',
 
 
 @entity_task(log, writes=['gridded_data'])
-def distribute_thickness_per_altitude(gdir, add_slope=True, topo='topo_smoothed',
+def distribute_thickness_per_altitude(gdir, add_slope=True,
+                                      topo_variable='topo_smoothed',
                                       smooth_radius=None,
                                       dis_from_border_exp=0.25,
                                       varname_suffix=''):
     """Compute a thickness map by redistributing mass along altitudinal bands.
 
-    This is a rather cosmetic task, not relevant for OGGM but for ITMIX.
+    This is a rather cosmetic task, not relevant for OGGM but for ITMIX or
+    for visualizations.
 
     Parameters
     ----------
     gdir : :py:class:`oggm.GlacierDirectory`
         the glacier directory to process
     add_slope : bool
         whether a corrective slope factor should be used or not
+    topo_variable : str
+        the topography to read from `gridded_data.nc` (could be smoothed, or
+        smoothed differently).
     smooth_radius : int
         pixel size of the gaussian smoothing. Default is to use
         cfg.PARAMS['smooth_window'] (i.e. a size in meters). Set to zero to
@@ -757,7 +762,7 @@ def distribute_thickness_per_altitude(gdir, add_slope=True, topo='topo_smoothed'
         gridded_attributes(gdir)
 
     with utils.ncDataset(grids_file) as nc:
-        topo_smoothed = nc.variables[topo][:]
+        topo = nc.variables[topo_variable][:]
         glacier_mask = nc.variables['glacier_mask'][:]
         dis_from_border = nc.variables['dis_from_border'][:]
         if add_slope:
@@ -786,25 +791,25 @@ def distribute_thickness_per_altitude(gdir, add_slope=True, topo='topo_smoothed'
 
     # Assign a first order thickness to the points
     # very inefficient inverse distance stuff
-    thick = glacier_mask * np.NaN
-    for y in range(thick.shape[0]):
-        for x in range(thick.shape[1]):
-            phgt = topo_smoothed[y, x]
-            # take the ones in a 100m range
-            starth = 100.
-            while True:
-                starth += 10
-                pok = np.nonzero(np.abs(phgt - hs) <= starth)[0]
-                if len(pok) != 0:
-                    break
-            sqr = np.sqrt((xs[pok]-x)**2 + (ys[pok]-y)**2)
-            pzero = np.where(sqr == 0)
-            if len(pzero[0]) == 0:
-                thick[y, x] = np.average(ts[pok], weights=1 / sqr)
-            elif len(pzero[0]) == 1:
-                thick[y, x] = ts[pzero]
-            else:
-                raise RuntimeError('We should not be there')
+    thick = glacier_mask * 0.
+    yglac, xglac = np.nonzero(glacier_mask == 1)
+    for y, x in zip(yglac, xglac):
+        phgt = topo[y, x]
+        # take the ones in a 100m range
+        starth = 100.
+        while True:
+            starth += 10
+            pok = np.nonzero(np.abs(phgt - hs) <= starth)[0]
+            if len(pok) != 0:
+                break
+        sqr = np.sqrt((xs[pok]-x)**2 + (ys[pok]-y)**2)
+        pzero = np.where(sqr == 0)
+        if len(pzero[0]) == 0:
+            thick[y, x] = np.average(ts[pok], weights=1 / sqr)
+        elif len(pzero[0]) == 1:
+            thick[y, x] = ts[pzero]
+        else:
+            raise RuntimeError('We should not be there')
 
     # Distance from border (normalized)
     dis_from_border = dis_from_border**dis_from_border_exp

oggm/core/massbalance.py

@@ -876,7 +876,8 @@ class RandomMassBalance(MassBalanceModel):
 
     def __init__(self, gdir, mb_model_class=MonthlyTIModel,
                  y0=None, halfsize=15, seed=None,
-                 all_years=False, unique_samples=False, prescribe_years=None,
+                 all_years=False, unique_samples=False,
+                 prescribe_years=None,
                  **kwargs):
         """Initialize.
 
@@ -1489,8 +1490,8 @@ def mb_calibration_from_geodetic_mb(gdir, *,
         ref_mb_df = get_geodetic_mb_dataframe().loc[gdir.rgi_id]
         ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == ref_period]
         # dmdtda: in meters water-equivalent per year -> we convert to kg m-2 yr-1
-        ref_mb = float(ref_mb_df['dmdtda']) * 1000
-        ref_mb_err = float(ref_mb_df['err_dmdtda']) * 1000
+        ref_mb = ref_mb_df['dmdtda'].iloc[0] * 1000
+        ref_mb_err = ref_mb_df['err_dmdtda'].iloc[0] * 1000
     except KeyError:
         if override_missing is None:
             raise