adding option for trapezoidal downstream line

oggm/cfg.py

@@ -526,6 +526,7 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
     PARAMS['mpi_recv_buf_size'] = cp.as_int('mpi_recv_buf_size')
     PARAMS['use_multiple_flowlines'] = cp.as_bool('use_multiple_flowlines')
     PARAMS['filter_min_slope'] = cp.as_bool('filter_min_slope')
+    PARAMS['downstream_line_shape'] = cp['downstream_line_shape']
     PARAMS['auto_skip_task'] = cp.as_bool('auto_skip_task')
     PARAMS['correct_for_neg_flux'] = cp.as_bool('correct_for_neg_flux')
     PARAMS['filter_for_neg_flux'] = cp.as_bool('filter_for_neg_flux')
@@ -595,7 +596,7 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
     except ValueError:
         PARAMS['prcp_scaling_factor'] = None
 
-     # Delete non-floats
+    # Delete non-floats
     ltr = ['working_dir', 'dem_file', 'climate_file', 'use_tar_shapefiles',
            'grid_dx_method', 'run_mb_calibration', 'compress_climate_netcdf',
            'mp_processes', 'use_multiprocessing', 'climate_qc_months',
@@ -616,7 +617,8 @@ def initialize_minimal(file=None, logging_level='INFO', params=None,
            'tidewater_type', 'store_model_geometry', 'use_winter_prcp_factor',
            'store_diagnostic_variables', 'store_fl_diagnostic_variables',
            'geodetic_mb_period', 'store_fl_diagnostics', 'winter_prcp_factor_ab',
-           'winter_prcp_factor_range', 'prcp_scaling_factor']
+           'winter_prcp_factor_range', 'prcp_scaling_factor',
+           'downstream_line_shape']
     for k in ltr:
         cp.pop(k, None)
 

oggm/core/centerlines.py

@@ -1103,10 +1103,12 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator):
         donot_compute.append(_isborder)
 
     bed = []
+    terrain_heights = []
     for ic, (cc, dontcomp) in enumerate(zip(cs, donot_compute)):
 
         if dontcomp:
             bed.append(np.NaN)
+            terrain_heights.append(np.NaN)
             continue
 
         z = []
@@ -1128,6 +1130,7 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator):
         extr = scipy.signal.argrelextrema(dsts, np.less, mode='wrap')
         if len(extr[0]) == 0:
             bed.append(np.NaN)
+            terrain_heights.append(np.NaN)
             continue
 
         # from local minima find that with the minimum |x|
@@ -1144,6 +1147,10 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator):
 
         p = _approx_parabola(roNx, zN, y0=zHead)
 
+        # define terrain height as the maximum height difference of points used
+        # for the parabolic fitting and the bottom height
+        terrain_heights.append(float(np.max(zN - zHead)))
+
         # shift parabola to the ds-line
         p2 = np.copy(p)
         p2[2] = z[ro == 0]
@@ -1157,6 +1164,9 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator):
         else:
             bed.append(np.NaN)
 
+    terrain_heights = np.asarray(terrain_heights)
+    assert len(terrain_heights) == idl.nx, 'len(terrain_heights) == idl.nx'
+
     bed = np.asarray(bed)
     assert len(bed) == idl.nx, 'len(bed) == idl.nx'
     pvalid = np.sum(np.isfinite(bed)) / len(bed) * 100
@@ -1169,6 +1179,8 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator):
     # interpolation, filling the gaps
     default = cfg.PARAMS['default_parabolic_bedshape']
     bed_int = interp_nans(bed, default=default)
+    default = 100.  # assume a default terrain height of 100 m if all NaN
+    terrain_heights = interp_nans(terrain_heights, default=default)
 
     # We forbid super small shapes (important! This can lead to huge volumes)
     # Sometimes the parabola fits in flat areas are very good, implying very
@@ -1178,12 +1190,42 @@ def _parabolic_bed_from_topo(gdir, idl, interpolator):
     # Smoothing
     bed_ma = pd.Series(bed_int)
     bed_ma = bed_ma.rolling(window=5, center=True, min_periods=1).mean()
-    return bed_ma.values
+    return bed_ma.values, terrain_heights
+
+
+def _trapezoidal_bottom_width_from_terrain_cross_section_area(
+        terrain_heights, Ps, lambdas, w0_min):
+    """This function calculates a bottom width for a trapezoidal downstream
+    line in a way that the terrain cross section area is preserved. The area to
+    preserve is defined by the fitted parabolic shape and the terrain height.
+
+    Parabolic formulas involved (Ap = area, h = terrain_height,
+    w = surface width, Ps = bed-shape parameter):
+        Ap = 2 / 3 * h * w
+        Ps = 4 * h / w^2
+
+    Trapezoidal formulas involved (At = area, h = terrain_height,
+    w = surface width, w0 = bottom width, lambda = defines angle of wall):
+        At = (w + w0) * h / 2
+        w = w0 + lambda * h
+
+    Putting all formulas together and setting Ap = At we get:
+        w0 = 4 / 3 * sqrt(h / Ps) - 1 / 2 * lambda * h
+    """
+
+    w0s = utils.clip_min(4 / 3 * np.sqrt(terrain_heights / Ps) -
+                         1 / 2 * lambdas * terrain_heights,
+                         w0_min)
+
+    return w0s
 
 
 @entity_task(log, writes=['downstream_line'])
 def compute_downstream_bedshape(gdir):
     """The bedshape obtained by fitting a parabola to the line's normals.
+    Further a trapezoidal shape is fitted to match the cross section area of
+    the valley. Which downstream shape (parabola or trapezoidal) is used can be
+    selected with cfg.PARAMS['downstream_line_shape'].
 
     Also computes the downstream's altitude.
 
@@ -1210,7 +1252,7 @@ def compute_downstream_bedshape(gdir):
     xy = (np.arange(0, len(y)-0.1, 1), np.arange(0, len(x)-0.1, 1))
     interpolator = RegularGridInterpolator(xy, topo)
 
-    bs = _parabolic_bed_from_topo(gdir, cl, interpolator)
+    bs, terrain_heights = _parabolic_bed_from_topo(gdir, cl, interpolator)
     assert len(bs) == cl.nx, 'len(bs) == cl.nx'
     assert np.all(np.isfinite(bs)), 'np.all(np.isfinite(bs))'
 
@@ -1222,10 +1264,20 @@ def compute_downstream_bedshape(gdir):
     # If smoothing, this is the moment
     hgts = gaussian_filter1d(hgts, cfg.PARAMS['flowline_height_smooth'])
 
+    # calculate bottom width of trapezoidal shapes
+    lambdas = np.ones(len(bs)) * cfg.PARAMS['trapezoid_lambdas']
+    w0_min = cfg.PARAMS['trapezoid_min_bottom_width']
+    w0s = _trapezoidal_bottom_width_from_terrain_cross_section_area(
+        terrain_heights, bs, lambdas, w0_min)
+    assert len(w0s) == cl.nx, 'len(w0s) == cl.nx'
+    assert np.all(np.isfinite(w0s)), 'np.all(np.isfinite(w0s))'
+    assert np.all(w0s >= w0_min), 'np.all(w0s >= w0_min)'
+
     # write output
     out = gdir.read_pickle('downstream_line')
     out['bedshapes'] = bs
     out['surface_h'] = hgts
+    out['w0s'] = w0s
     gdir.write_pickle(out, 'downstream_line')
 
 

oggm/core/flowline.py

@@ -2787,14 +2787,26 @@ def init_present_time_glacier(gdir, filesuffix=''):
         lambdas[bed_shape == 0] = def_lambda
 
         if not gdir.is_tidewater and inv['is_last']:
-            # for valley glaciers, simply add the downstream line
+            # for valley glaciers, simply add the downstream line, depending on
+            # selected shape parabola or trapezoidal
             dic_ds = gdir.read_pickle('downstream_line')
-            bed_shape = np.append(bed_shape, dic_ds['bedshapes'])
-            lambdas = np.append(lambdas, dic_ds['bedshapes'] * np.NaN)
+            if cfg.PARAMS['downstream_line_shape'] == 'parabola':
+                bed_shape = np.append(bed_shape, dic_ds['bedshapes'])
+                lambdas = np.append(lambdas, dic_ds['bedshapes'] * np.NaN)
+                widths_m = np.append(widths_m, dic_ds['bedshapes'] * 0.)
+            elif cfg.PARAMS['downstream_line_shape'] == 'trapezoidal':
+                bed_shape = np.append(bed_shape, dic_ds['bedshapes'] * np.NaN)
+                lambdas = np.append(lambdas, np.ones(len(dic_ds['w0s'])) *
+                                    def_lambda)
+                widths_m = np.append(widths_m, dic_ds['w0s'])
+            else:
+                raise InvalidParamsError(
+                    f"Unknown cfg.PARAMS['downstream_line_shape'] = "
+                    f"{cfg.PARAMS['downstream_line_shape']} (options are "
+                    f"'parabola' and 'trapezoidal').")
             section = np.append(section, dic_ds['bedshapes'] * 0.)
             surface_h = np.append(surface_h, dic_ds['surface_h'])
             bed_h = np.append(bed_h, dic_ds['surface_h'])
-            widths_m = np.append(widths_m, dic_ds['bedshapes'] * 0.)
             line = dic_ds['full_line']
 
         if gdir.is_tidewater and inv['is_last']:

oggm/params.cfg

@@ -164,6 +164,14 @@ base_binsize = 50.
 # 1 means default (i.e. kernel size 9).
 smooth_widths_window_size = 1
 
+### DOWNSTREAM LINE
+# define which bed shape should be used for the ice free downstream line, the
+# options are 'parabola' or 'trapezoidal'
+downstream_line_shape = 'parabola'
+# defines the minimum bottom width of a potential trapezoidal downstream line
+# in meters
+trapezoid_min_bottom_width = 50.
+
 ### CLIMATE params
 
 # Baseline climate is the reference climate data to use for this workflow.

oggm/tests/test_models.py

@@ -67,9 +67,11 @@
 
 class TestInitPresentDayFlowline:
 
-    def test_init_present_time_glacier(self, hef_gdir):
+    @pytest.mark.parametrize('downstream_line_shape', ['parabola', 'trapezoidal'])
+    def test_init_present_time_glacier(self, hef_gdir, downstream_line_shape):
 
         gdir = hef_gdir
+        cfg.PARAMS['downstream_line_shape'] = downstream_line_shape
         init_present_time_glacier(gdir)
 
         fls = gdir.read_pickle('model_flowlines')
@@ -116,6 +118,19 @@ def test_init_present_time_glacier(self, hef_gdir):
         np.testing.assert_allclose(gdir.rgi_area_km2, ref_area * 1e-6)
         np.testing.assert_allclose(gdir.rgi_area_km2, area)
 
+        # test that final downstream line has the desired shape
+        fl = fls[-1]
+        ice_mask = fls[-1].thick > 0.
+        if downstream_line_shape == 'parabola':
+            assert np.all(np.isfinite(fl.bed_shape[~ice_mask]))
+            assert np.all(~np.isfinite(fl._w0_m[~ice_mask]))
+        if downstream_line_shape == 'trapezoidal':
+            assert np.all(np.isfinite(fl._w0_m[~ice_mask]))
+            assert np.all(~np.isfinite(fl.bed_shape[~ice_mask]))
+            # check that bottom width of downstream line is larger than minimum
+            assert np.all(fl._w0_m[~ice_mask] >
+                          cfg.PARAMS['trapezoid_min_bottom_width'])
+
         if do_plot:
             plt.plot(fls[-1].bed_h, color='k')
             plt.plot(fls[-1].surface_h)
@@ -127,6 +142,10 @@ def test_init_present_time_glacier(self, hef_gdir):
         init_present_time_glacier(gdir, filesuffix='_test')
         assert os.path.isfile(os.path.join(gdir.dir, 'model_flowlines_test.pkl'))
 
+        cfg.PARAMS['downstream_line_shape'] = 'free_shape'
+        with pytest.raises(InvalidParamsError):
+            init_present_time_glacier(gdir)
+
     def test_present_time_glacier_massbalance(self, hef_gdir):
 
         gdir = hef_gdir