Non-linear bed slopes and collection updates. (#139)

* Passing partial expressions as values in change_attributes dict.

This makes is possible to pass partial mathematical expressions when
updating attributes through the change_attributes method of
GlacierCollection. This is done by passing a string e.g. '* 10' in the
value list of a key, This will multiply the current value by a factor of
10.

Partial expressions are evaluated in the new utility function expression_parser.

Adds test cases for both the expression_parser and the change_attributes
method.

* Non-linear bed profiles.

This give the user more control over the slope of the glacier bed.
Either by supplying only a single slope or a sequence of slopes and
slope sections (specified by altitude intervals).

Behind the scenes it figures out how long (x-distance) each section has
to be to fulfill the slope, which is what the total glacier length and
the number of gridpoints is then based on. It then interpolates the
slope sections along the distance of the glacier.

* First rules to _check_collection.

With the introduction of non-linear bed slopes, the collection need
logic to determine if glaciers can be put in the same collection. In
this initial case we check if the bed profiles are the same. This is to
keep plot methods working. Maybe this is something will do different in
the future. But for now I don't see a nice way to plot glaciers with
different slopes on top of each other.

* Added checks to make sure provided slopes are reasonable.

oggm_edu/funcs.py

@@ -79,11 +79,54 @@ def context_wrapper(
         figsize=(12, 9),
         **kwargs,
     ):
-        with \
-                mpl.rc_context({"figure.figsize": figsize}), \
-                sns.plotting_context(sns_context), \
-                sns.axes_style(sns_axes_style)\
-                :
+        with mpl.rc_context({"figure.figsize": figsize}), sns.plotting_context(
+            sns_context
+        ), sns.axes_style(sns_axes_style):
             return func(*args, **kwargs)
 
     return context_wrapper
+
+
+def expression_parser(expression, numeric):
+    """Parse an uncomplete math expression e.g. '* 10' and apply it to supplied attribute.
+
+    Parameters
+    ----------
+    expression : string
+        Incomplete math expression in the form '* 5' or '- 10'. Can also be empty to leave
+        numeric un-affected.
+    numeric : int or float
+        Value to evaluate against expression.
+    Returns
+    -------
+    The full expression evaluated.
+    """
+
+    # is expression a string?
+    if not isinstance(expression, str):
+        raise TypeError("expression should be a string.")
+
+    elif not isinstance(numeric, (int, float)):
+        raise TypeError("numeric should be an integer or a float.")
+    # If expression is empty, we return the numeric.
+    elif expression == "":
+        return numeric
+
+    else:
+        # Extract the operator
+        operator = expression[0]
+        if operator not in ["*", "/", "+", "-"]:
+            raise ValueError(
+                "First part of expression should be either one of *, /, +, -."
+            )
+        # Extract the factor
+        factor = float(expression[1:])
+        # Create a table of possible expressions. Then we just pick the correct one for return.
+        expression_dict = {
+            "*": numeric * factor,
+            "/": numeric / factor,
+            "+": numeric + factor,
+            "-": numeric - factor,
+        }
+
+        return expression_dict[operator]

oggm_edu/glacier_bed.py

@@ -3,6 +3,7 @@
 """
 
 from oggm_edu.funcs import edu_plotter
+from collections.abc import Sequence
 
 # Other libraries.
 import numpy as np
@@ -22,14 +23,18 @@ def __init__(
         width=None,
         altitudes=None,
         widths=None,
-        slope=None,
+        slopes=None,
+        slope_sections=None,
         nx=200,
         map_dx=100,
     ):
         """Initialise the bed. Pass single scalars for top, bottom and width
-         to create a square glacier bed. For finer control pass altitudes and
-         widths as lists or tuples for custom geometry. Will linearly
-         intepolate between the altitude/width pairs.
+        to create a square glacier bed. For finer control of the width pass
+        altitudes and widths as lists or tuples for custom geometry. This will linearly
+        intepolate between the altitude/width pairs.
+
+        To contol the slope of the glacier provide a single value to slopes. For even finer contol
+        pass a sequence of values to slopes along with a sequence of altitudes to slope_sections.
 
         Parameters
         ----------
@@ -46,10 +51,15 @@ def __init__(
         widths : array_like, (ints of floats)
             List of values defining the widths along the glacier.
             Length should match altitudes.
-        slope : float
-            Define the slope of the bed, decimal
+        slopes : array_like(int)
+            Define the slope of the bed, degrees. One or multiple values. For the latter,
+            slope sections are required.
+        slope_sections : array_like(int)
+            Defines the altitude spans for the slopes. Should start with the top altitude and
+            end with the bottom altitude. The number of altitudes need to be one greater then
+            the number of slopes.
         nx : int
-            Number of grid points.
+            Number of grid points. Will be overriden if slope is provided.
         map_dx : int
             Grid point spacing in meters.
         """
@@ -95,30 +105,80 @@ def __init__(
                 "Top of the bed has to be above the bottom."
                 + " Bottom also has to be above 0"
             )
-        # Calculate the slope correction
-        if slope:
-            if slope >= 1.0 or slope <= 0.0:
-                raise ValueError("Slope shoud be between 0 and 1 (not equal)")
-            # If there is a slop we re calculate the map_dx
-            else:
-                map_dx = (self.top - self.bottom) / (nx * slope)
-
         # Set the resolution.
         self.map_dx = map_dx
-        # Initialise the bed
-        self.bed_h = np.linspace(self.top, self.bottom, nx)
-        self.distance_along_glacier = np.linspace(0, nx, nx) * self.map_dx * 1e-3
+        # Do we have a specified slope?
+        if slopes:
+            # If slopes is not a sequence, then make it one.
+            if not isinstance(slopes, Sequence):
+                slopes = [slopes]
+
+            # Make sure that the provided slopes are reasonable.
+            if (np.asarray(slopes) < 1).any() or (np.asarray(slopes) > 80).any():
+                raise ValueError("Slopes should be above 0 and below 80 degrees.")
+            #  Do we have sequence of both slopes and breakpoints?
+            if isinstance(slopes, Sequence) and isinstance(slope_sections, Sequence):
+                # Are they compatible?
+                # There should be one more section compared to slopes
+                if not len(slopes) == len(slope_sections) - 1:
+                    raise ValueError(
+                        "Number of slope sections should be one more then number of slopes"
+                    )
+                # Have to match top and bottom.
+                elif slope_sections[0] != self.top or slope_sections[-1] != self.bottom:
+                    raise ValueError(
+                        "First and last value of slope_sections should match top and bottom."
+                    )
+            # If we passed a single slope, we can assign slope sections to still make use of our fancy algo.
+            elif len(slopes) == 1:
+                slope_sections = [self.top, self.bottom]
+            # What is the height difference between the sections?
+            slope_sections = np.asarray(slope_sections)
+            slope_sections_diff = np.abs(np.diff(slope_sections))
+
+            # How long does a segment has to be to have the correct slope angle?
+            x_segments_length = slope_sections_diff / np.tan(np.deg2rad(slopes))
+            # We add a zero to the start to begin the interpolation in the right place.
+            # Take the cumsum to get the absolute position of each coord.
+            x_segments_length = np.concatenate([[0], x_segments_length.cumsum()])
+
+            # What does the total length of the glacier have to be?
+            # This should be a float for correct interpolation.
+            total_length = x_segments_length.max()
+            # And how many gridpoints do we need for this?
+            self.nx = int(total_length / self.map_dx)
+            # Distance along the glacier: up to the total length in the closest number of grid points.
+            self.distance_along_glacier = np.linspace(0, total_length, self.nx)
+            # Interpolate the heights
+            heights = np.interp(
+                self.distance_along_glacier, x_segments_length, slope_sections
+            )
+
+            # We can now put the distance in kms
+            self.distance_along_glacier = self.distance_along_glacier * 1e-3
+
+            # Assign the heights.
+            self.bed_h = heights
+
+        # Otherwise we just make a simple bed profile.
+        else:
+            self.nx = nx
+            self.bed_h = np.linspace(self.top, self.bottom, self.nx)
+            self.distance_along_glacier = (
+                np.linspace(0, self.nx, self.nx) * map_dx * 1e-3
+            )
+
         # Widths.
         # If width has a length of one, we have a constant width.
         if width:
             self.width = width
-            self.widths = (np.zeros(nx) + self.width) / map_dx
+            self.widths = (np.zeros(self.nx) + self.width) / self.map_dx
         # If length of width and length of width altitudes are compatible,
         # we create a bed with variable width.
         elif len(altitudes) == len(widths):
             self.width = widths
             # First create a constant bed.
-            tmp_w = np.zeros(nx)
+            tmp_w = np.zeros(self.nx)
 
             # Make sure we have lists.
             altitudes = list(altitudes)
@@ -138,7 +198,7 @@ def __init__(
                 # Linear interpolation between the widths.
                 tmp_w[mask] = np.linspace(self.width[i], self.width[i + 1], mask.sum())
             # Assign the varied widths it.
-            self.widths = tmp_w / map_dx
+            self.widths = tmp_w / self.map_dx
 
         # If noting works, raise an error.
         else:

oggm_edu/glacier_collection.py

@@ -4,7 +4,7 @@
 
 # Internals
 from oggm_edu.glacier import Glacier
-from oggm_edu.funcs import edu_plotter
+from oggm_edu.funcs import edu_plotter, expression_parser
 
 # Other libraries.
 import pandas as pd
@@ -27,15 +27,15 @@ def __init__(self, glacier_list=None):
         ----------
         glacier_list : list of glaciers objects
             Defaults to none. Has the possibility to add glaciers on the go.
+            Glaciers need to have the same bed slope.
         """
 
         self._glaciers = []
 
         # Do we have a list of glaciers?
         if glacier_list:
-            # Check that all are Glaciers.
-            if all(isinstance(glacier, Glacier) for glacier in glacier_list):
-                self._glaciers = glacier_list
+            # Use the add method.
+            self.add(glacier_list)
 
     def __repr__(self):
         return f"Glacier collection with {len(self._glaciers)} glaciers."
@@ -70,10 +70,21 @@ def _repr_html_(self):
         else:
             pass
 
-    def check_collection(self, glacier):
+    def _check_collection(self, glacier):
         """Utility method. Check if the glaciers obey the collection rules.
-        They need to have the same domains for this to make sense I think."""
-        # TODO
+        Make sure that a new glacier has the same bed as other glaciers in the collection.
+
+        glacier : Glacier
+            New glacier to check against the collection.
+        """
+
+        # Since this happens as soon as a glacier is added to a collection, we only have to check
+        # against the previous glacier
+        beds_ok = np.array_equal(self._glaciers[-1].bed.bed_h, glacier.bed.bed_h)
+
+        # If beds are not ok
+        if not beds_ok:
+            raise ValueError("Bed of new glacier does not match beds in collection.")
 
     @property
     def glaciers(self):
@@ -116,7 +127,7 @@ def fill(self, glacier, n, attributes_to_change=None):
             self.change_attributes(attributes_to_change)
 
     def add(self, glacier):
-        """Add one or more glaciers to the collection.
+        """Add one or more glaciers to the collection. Glaciers have to have the same slope.
 
         Parameters
         ----------
@@ -132,24 +143,26 @@ def add(self, glacier):
                         "Glacier collection can only contain glacier objects."
                     )
                 # Is the glacier already in the collection?
-                if glacier in self._glaciers:
+                elif glacier in self._glaciers:
                     raise AttributeError("Glacier is already in collection")
-                # If the glacier is an instance of Glacier, we can add it to
-                # the collection.
-                else:
-                    self._glaciers.append(glacier)
+                # If there are already glaciers in the collection, check that the new one fit.
+                elif self._glaciers:
+                    self._check_collection(glacier)
+                # If no throws, add it.
+                self._glaciers.append(glacier)
         # If not iterable
         else:
             # Check that glacier is of the right type.
             if not isinstance(glacier, Glacier):
                 raise TypeError("Glacier collection can only contain glacier objects.")
             # Is the glacier already in the collection?
-            if glacier in self._glaciers:
+            elif glacier in self._glaciers:
                 raise AttributeError("Glacier is already in collection")
-            # If the glacier is an instance of Glacier, we can add it to
-            # the collection.
-            else:
-                self._glaciers.append(glacier)
+            # If there are already glaciers in the collection, check that the new one fit.
+            elif self._glaciers:
+                self._check_collection(glacier)
+            # If no throws, add it.
+            self._glaciers.append(glacier)
 
     def change_attributes(self, attributes_to_change):
         """Change the attribute(s) of the glaciers in the collection.
@@ -160,14 +173,18 @@ def change_attributes(self, attributes_to_change):
             Dictionary where the key value pairs follow the structure:
             {'key': [n values], ...}, where 'key' match an attribute
             of the glacier and n match the length of the collection.
-            Valid values for key are:
+            Valid keys are:
             - gradient
             - ela
             - basal_sliding
             - creep
             - normal_years
             - surging_years
             - basal_sliding_surge
+            Values should be either numeric or a string of a partial
+            mathematical expression e.g. '* 10', which is evaluated
+            to multiplying the current value by a factor 10. Pass an empty
+            string to leave the attribute unaffected.
         """
         # Did we get a dict?
         if not isinstance(attributes_to_change, dict):
@@ -198,13 +215,26 @@ def change_attributes(self, attributes_to_change):
                 )
             # If all passes
             else:
+                # Set values and glaciers.
                 for (glacier, value) in zip(self.glaciers, values):
-                    # Use built in setattr. Should respect the defined
-                    # setters, with error messages an such.
+                    # Should we act on the glacier or mass balance?
                     if key in mb_attrs:
-                        setattr(glacier.mass_balance, key, value)
+                        obj = glacier.mass_balance
+                    # Just the glacier
                     else:
-                        setattr(glacier, key, value)
+                        obj = glacier
+
+                    # If value is a string (a partial expression) we set the value
+                    # using the expression_parser.
+                    if isinstance(value, str):
+                        # Get the current value of the attribute.
+                        curr_value = getattr(obj, key)
+                        # Calculate the value based on the partial expression.
+                        value = expression_parser(value, curr_value)
+                    # If value is not a string we can just assign it directly.
+                    # Use built in setattr. Should respect the defined
+                    # setters, with error messages an such.
+                    setattr(obj, key, value)
 
     def progress_to_year(self, year):
         """Progress the glaciers within the collection to

oggm_edu/tests/test_funcs.py

@@ -1,6 +1,8 @@
 import oggm_edu
+from oggm_edu.funcs import expression_parser
 import matplotlib.pyplot as plt
 from numpy.testing import assert_equal
+import pytest
 
 
 def test_plot_glacier_graphics():
@@ -18,3 +20,19 @@ def plot_sth():
 
     ax = plot_sth(figsize=(16, 16))
     assert_equal(ax.get_figure().get_size_inches(), [16.0, 16.0])
+
+
+def test_expression_parser():
+    """Test the expression parser."""
+    # Does it work like we think it does?
+    assert 20 == expression_parser("*10", 2)
+    assert 4 == expression_parser("+ 2", 2.0)
+    assert 4.5 == expression_parser("+ 2.5", 2.0)
+    assert 1 == expression_parser("/ 2", 2.0)
+    assert 5 == expression_parser("", 5.0)
+
+    # Should raise
+    with pytest.raises(Exception) as e_info:
+        expression_parser(" * 10", 2)
+    with pytest.raises(Exception) as e_info:
+        expression_parser(" * 10", "elk")