Add documentation for inversion-task

docs/_code/apply_invert_1.py

@@ -0,0 +1,2 @@
+from oggmcontrib.tasks import distributed_vas_thickness
+out_thick = distributed_vas_thickness(gdir)

docs/_code/apply_invert_2.py

@@ -0,0 +1,2 @@
+from oggmcontrib.tasks import distributed_vas_thickness
+out_thick = distributed_vas_thickness(gdir, sqrt=True)

docs/_code/invert_multi_plotfunc.py

@@ -0,0 +1,11 @@
+import netCDF4
+
+def plot_inversion(gdir, ax):
+    """Plot the VAS inversion for this glacier."""
+
+    # Read the data
+    grids_file = gdir.get_filepath('gridded_data')
+    with netCDF4.Dataset(grids_file) as nc:
+        thick = nc.variables['vas_distributed_thickness'][:]
+
+    ax.imshow(thick)

docs/_code/prepro_invert.py

@@ -0,0 +1,18 @@
+import geopandas as gpd
+import oggm
+from oggm import cfg, tasks
+from oggm.utils import get_demo_file
+
+# Set up the input data for this example
+cfg.initialize()
+cfg.PATHS['working_dir'] = oggm.gettempdir('oggm_wd')
+cfg.PATHS['dem_file'] = get_demo_file('srtm_oetztal.tif')
+cfg.set_intersects_db(get_demo_file('rgi_intersect_oetztal.shp'))
+
+# Glacier directory for Hintereisferner in Austria
+entity = gpd.read_file(get_demo_file('Hintereisferner_RGI5.shp')).iloc[0]
+gdir = oggm.GlacierDirectory(entity)
+
+# The usual OGGM preprecessing
+tasks.define_glacier_region(gdir, entity=entity)
+tasks.glacier_masks(gdir)

docs/_code/prepro_invert_multi.py

@@ -0,0 +1,10 @@
+# Read in the RGI file
+rgi_file = get_demo_file('rgi_oetztal.shp')
+rgidf = gpd.read_file(rgi_file)
+
+# Use multiprocessing to apply the OGGM tasks and the new task to all glaciers
+from oggm import workflow
+gdirs = workflow.init_glacier_regions(rgidf)
+workflow.execute_entity_task(tasks.glacier_masks, gdirs)
+# Yes, also your new task!
+workflow.execute_entity_task(distributed_vas_thickness, gdirs)

docs/index.rst

@@ -35,22 +35,28 @@ How does this package work?
 
 This package is a template. It implements two very simple use cases:
 
-- adding a new "bed inversion" task based on volume-area-scaling but using the
-  OGGM preprocessing workflow
-- using a custom mass-balance model and apply it to OGGM's ice-dynamics model
+- :ref:`inversion-task` based on volume-area-scaling
+  but using the OGGM preprocessing workflow
+- :ref:`mass-balance` and apply it to OGGM's ice dynamics model
 
-You can install this custom package with:
+
+Installation, usage
+~~~~~~~~~~~~~~~~~~~
+
+You can install this custom package with::
 
      pip install git+https://github.com/OGGM/oggmcontrib.git
 
-However, what you'd probably like to do is to `fork <https://help.github.com/articles/fork-a-repo/>`_ this repository and use
-it as a template for your own project. You can install it locally with
+However, what you'd probably like to do is to `fork <https://help.github.com/articles/fork-a-repo/>`_
+this repository and use it as a template for your own project.
+
+Once available locally, you can install it in development mode with::
 
     pip install -e .
 
 
-Examples
-~~~~~~~~
+Index
+~~~~~
 
 .. toctree::
     :maxdepth: 1

docs/inversion-task.rst

@@ -1,3 +1,140 @@
+.. _inversion-task:
+
 Adding an ice thickness inversion task
 ======================================
 
+This example illustrates the concept of `entity tasks`_ by implementing
+a new task for the ice thickness inversion. We are therefore using
+the OGGM preprocessing workflow only, and add our own tool on top of it.
+
+The new task is called ``distributed_vas_thickness`` and is very simple:
+with a baseline volume obtained from volume area scaling, we distribute the
+local ice thickness so that the ice gets thicker as a function of the
+distance to the boundaries, and so that the total volume is conserved.
+
+The implementation can be found in `oggmcontrib/tasks.py`_. This module contains
+two functions: a helper function called ``distance_from_border`` and the
+actual task ``distributed_vas_thickness``.
+
+Have a short look at the code in `oggmcontrib/tasks.py`_ before going on.
+
+**Entity tasks** in OGGM allow only one single argument: a `GlacierDirectory`_
+instance, which gives access to all the input files in the working directory.
+Tasks can have as many keyword arguments as you wish, though. They can return
+data (useful for testing), but in order to fit in the stop/restart workflow
+of OGGM they should write their output in the working directory though. Here
+we use an existing NetCDF file and add the output of our function to it.
+
+The last element that makes of a function a real "entity task" for OGGM is the
+addition of the ``@entity_task`` decorator on top of it. If you are not used
+to python decorators, don't worry: just keep in mind that these decorators are
+here for three major purposes:
+
+- logging
+- error handling (if your task raises an error on certain glaciers, OGGM might
+  choose to ignore it if the user wants it that way)
+- multi-processing (by complying to a certain syntax, tasks cn be sent to
+  OGGM's task manager)
+
+.. _entity tasks: http://oggm.readthedocs.io/en/latest/api.html#entity-tasks
+.. _oggmcontrib/tasks.py: https://github.com/OGGM/oggmcontrib/blob/master/oggmcontrib/tasks.py
+.. _GlacierDirectory: http://oggm.readthedocs.io/en/latest/generated/oggm.GlacierDirectory.html#oggm.GlacierDirectory
+
+
+Apply our new task to a single glacier
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. ipython:: python
+   :suppress:
+
+   fpath = "_code/prepro_invert.py"
+   with open(fpath) as f:
+      code = compile(f.read(), fpath, 'exec')
+      exec(code)
+
+Let's use our usual test glacier for this:
+
+.. literalinclude:: _code/prepro_invert.py
+
+This was all OGGM stuff. Now let's use our new task on this preprocessed data:
+
+.. ipython:: python
+   :suppress:
+
+   fpath = "_code/apply_invert_1.py"
+   with open(fpath) as f:
+      code = compile(f.read(), fpath, 'exec')
+      exec(code)
+
+.. literalinclude:: _code/apply_invert_1.py
+
+
+And plot it:
+
+.. ipython:: python
+
+   plt.figure(figsize=(7, 4));
+   plt.imshow(out_thick);
+   @savefig plot_thick_1.png width=80%
+   plt.colorbar(label='Thick [m]');
+
+
+We can use the keyword arguments just like a regular function of course:
+
+.. ipython:: python
+   :suppress:
+
+   fpath = "_code/apply_invert_2.py"
+   with open(fpath) as f:
+      code = compile(f.read(), fpath, 'exec')
+      exec(code)
+
+.. literalinclude:: _code/apply_invert_2.py
+
+.. ipython:: python
+
+   plt.figure(figsize=(7, 4));
+   plt.imshow(out_thick);
+   @savefig plot_thick_2.png width=80%
+   plt.colorbar(label='Thick [m]');
+
+
+Apply our new task in parallel
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. ipython:: python
+   :suppress:
+
+   fpath = "_code/prepro_invert_multi.py"
+   with open(fpath) as f:
+      code = compile(f.read(), fpath, 'exec')
+      exec(code)
+
+Let's go big, and apply our task to the selections of glaciers in the
+Öztal Alps:
+
+.. literalinclude:: _code/prepro_invert_multi.py
+
+Define a simple function to plot them:
+
+.. ipython:: python
+   :suppress:
+
+   fpath = "_code/invert_multi_plotfunc.py"
+   with open(fpath) as f:
+      code = compile(f.read(), fpath, 'exec')
+      exec(code)
+
+.. literalinclude:: _code/invert_multi_plotfunc.py
+
+
+Let's go:
+
+.. ipython:: python
+
+
+   f, axs = plt.subplots(3, 3, figsize=(7, 7));
+   for gdir, ax in zip(gdirs[:9], np.array(axs).flatten()):
+      plot_inversion(gdir, ax)
+   @savefig plot_thick_all.png width=100%
+   plt.tight_layout();

docs/mass-balance.rst

@@ -1,3 +1,5 @@
+.. _mass-balance:
+
 Adding a mass-balance model
 ===========================
 

oggmcontrib/tasks.py

@@ -2,6 +2,7 @@
 
 """
 import logging
+import warnings
 import numpy as np
 import netCDF4
 import pandas as pd
@@ -17,11 +18,22 @@
 log = logging.getLogger(__name__)
 
 
-def distance_from_border(gdir):
+def distance_from_border(gdir,  sqrt=False):
     """Computes a normalized distance from border mask.
 
     It is a bit cleverer than just taking the glacier outlines: we also
     consider where the glacier has neighbors.
+
+    Parameters
+    ----------
+    gdir : oggm.GlacierDirectory
+        the working glacier directory
+    sqrt : bool, optional
+        whether the square root of the distance mask should be used or not
+
+    Returns
+    -------
+    the distance
     """
 
     # Variables
@@ -50,7 +62,9 @@ def distance_from_border(gdir):
     jj, ii = np.where(glacier_ext)
     for j, i in zip(jj, ii):
         dist = np.append(dist, np.min(gdfi.distance(shpg.Point(i, j))))
-    pok = np.where(dist <= 1)
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+        pok = np.where(dist <= 1)
     glacier_ext_intersect = glacier_ext * 0
     glacier_ext_intersect[jj[pok], ii[pok]] = 1
 
@@ -59,35 +73,47 @@ def distance_from_border(gdir):
     dis_from_border = distance_transform_edt(dis_from_border) * dx
     dis_from_border[np.where(glacier_mask == 0)] = np.NaN
 
+    # Square?
+    if sqrt:
+        dis_from_border = np.sqrt(dis_from_border)
+
     # We normalize and return
-    return dis_from_border / dis_from_border.sum()
+    return dis_from_border / np.nansum(dis_from_border)
 
 
 # An entity task decorator is what allows this task to work like other
 # OGGM tasks. What's happening inside is your job!
+# Here we document our task by saying that it will write in the gridded_data
+# file, but this is not required
 @entity_task(log, writes=['gridded_data'])
-def distributed_vas_thickness(gdir):
+def distributed_vas_thickness(gdir, sqrt=False):
     """Compute a thickness map of the glacier using very simple rules.
 
     The rules are:
     - the volume of the glacier is obtained from Volume Area Scaling (VAS)
     - the glacier gets thicker as a function of the distance to the outlines.
+
+    Parameters
+    ----------
+    gdir : oggm.GlacierDirectory
+        the working glacier directory
+    sqrt : bool, optional
+        whether the square root of the distance mask should be used or not
     """
 
     # Variables
-    dis_from_border = distance_from_border(gdir)
-    mask = np.isfinite(dis_from_border)
+    dis_from_border = distance_from_border(gdir,
+                                           sqrt=sqrt)
 
-    # VAS volume
+    # VAS volume in m3
     inv_vol = 0.034 * (gdir.rgi_area_km2**1.375)
+    inv_vol *= 1e9
 
     # Naive thickness
     dx = gdir.grid.dx
-    area_m2 = np.nansum(mask * dx**2)
-    thick = mask * inv_vol / area_m2
+    thick = dis_from_border * inv_vol
 
     # Conserve volume
-    # (by construction this shouldn't be necessary, but numerical errors...)
     thick *= inv_vol / np.nansum(thick * dx**2)
 
     # Write - we add it to existing netCDF file