Add profile interpolation method.

doc/api/datasets.rst

@@ -21,6 +21,7 @@ Postprocessing
    assign_isostasy
    getvar
    interp
+   profile
    where
    where_icemask
    where_thicker

doc/conf.py

@@ -60,6 +60,7 @@
 intersphinx_mapping = {
     'geopandas': ('https://geopandas.org/en/stable/', None),
     'matplotlib': ('https://matplotlib.org/stable/', None),
+    'numpy': ('https://numpy.org/doc/stable/', None),
     'pandas': ('https://pandas.pydata.org/docs/', None),
     'python': ('https://docs.python.org/3/', None),
     'xarray': ('http://xarray.pydata.org/en/stable/', None)

doc/datasets/masking.rst

@@ -4,8 +4,8 @@
 Masking and interpolation
 =========================
 
-Ice thickness threshold
------------------------
+Masking variables
+-----------------
 
 Hyoga's plot methods use an ice mask to determine which grid cells are
 glacierized and which are not. According to CF conventions, this is defined by
@@ -33,10 +33,6 @@ the case in the demo files.
    # restore the default of 1 m
    hyoga.config.glacier_masking_point = 1
 
-
-Custom glacier mask
--------------------
-
 For more control, on can set the ``land_ice_area_fraction`` variable using
 :meth:`~.Dataset.hyoga.assign_icemask`. Suppose that we define glaciers as grid
 cells filled with ice at least a metre thick, and moving at least ten metres
@@ -65,8 +61,8 @@ values with `np.nan` outside the where condition. However, they are meant to
 only affect "glacier variables" (currently any variable whose standard name
 does not start with ``bedrock_altitude``).
 
-Interpolation to another topography
------------------------------------
+Grid interpolation
+------------------
 
 For enhanced visuals, hyoga supports interpolating and remasking model output
 onto a different (usually higher-resolution) topography. This produces an image
@@ -119,5 +115,67 @@ with a much higher resolution.
    ds.hyoga.plot.ice_margin(edgecolor='0.25')
    ds.hyoga.plot.scale_bar()
 
-.. _xarray: https//xarray.pydata.org
-.. _`CF standard names`: http://cfconventions.org/standard-names.html
+Profile interpolation
+---------------------
+
+Profile interpolation aggregates variables with two horizontal dimensions (and
+possibly more dimensions) along a profile curve (or surface) defined by a
+sequence of (x, y) points. Let us draw a linear cross-section across the
+example dataset and plot this line on a map.
+
+.. plot::
+   :context: reset
+
+   # prepare a linear profile
+   import numpy as np
+   x = np.linspace(250e3, 450e3, 21)
+   y = np.linspace(5200e3, 5000e3, 21)
+   points = list(zip(x, y))
+
+   # plot profile line on a map
+   ds = hyoga.open.example('pism.alps.out.2d.nc')
+   ds.hyoga.plot.bedrock_altitude(center=False)
+   ds.hyoga.plot.ice_margin(facecolor='tab:blue')
+   plt.plot(x, y, color='tab:red', marker='x')
+
+The accessor method :meth:`~.Dataset.hyoga.profile` can be used to linearly
+interpolate the gridded dataset onto these coordinates, producing a new dataset
+where the ``x`` and ``y`` coordinates are swapped for a new coordinate ``d``,
+for the distance along the profile.
+
+.. plot::
+   :context: close-figs
+
+   profile = ds.hyoga.profile(points)
+   profile.hyoga.getvar('bedrock_altitude').plot(color='0.25')
+   profile.hyoga.getvar('surface_altitude').plot(color='C0')
+
+An additional ``interval`` keyword can be passed to control the horizontal
+resolution of the new profile dataset. This is particularly useful if the
+sequence of ``points`` is not regularly spaced.
+
+.. plot::
+   :context: close-figs
+
+   # prepare three subplots
+   fig, axes = plt.subplots(nrows=3, sharex=True, sharey=True)
+
+   # 10, 3, and 1 km resolutions
+   for ax, interval in zip(axes, [10e3, 3e3, 1e3]):
+       profile = ds.hyoga.profile(points, interval=interval)
+       profile.hyoga.getvar('bedrock_altitude').plot(ax=ax, color='0.25')
+       profile.hyoga.getvar('surface_altitude').plot(ax=ax, color='C0')
+
+   # remove duplicate labels
+   for ax in axes[:2]:
+       ax.set_xlabel('')
+   for ax in axes[::2]:
+       ax.set_ylabel('')
+
+The sequence of points in the above example does not have to form a straight
+line. Besides, it can also be provided as a :class:`numpy.ndarray`, a
+:class:`geopandas.GeoDataFrame`, or a path to a vector file containing a
+single, linestring geometry. Here is a more advanced example using a custom
+shapefile provided in the example data.
+
+.. plot:: ../examples/interp/plot_profile_altitudes.py

doc/roadmap.rst

@@ -56,7 +56,7 @@ v0.2.x Cartography
 Datasets
 ~~~~~~~~
 
-- |-| :meth:`.Dataset.hyoga.profile`
+- |x| :meth:`.Dataset.hyoga.profile`
 - |x| :meth:`.Dataset.hyoga.plot.bedrock_hillshade`
 - |x| :meth:`.Dataset.hyoga.plot.surface_hillshade`
 - |x| :meth:`.Dataset.hyoga.plot.natural_earth`
@@ -76,7 +76,6 @@ Documentation
 - |-| doc: ``foreword/history``
 - |x| :doc:`datasets/shading`
 - |x| :doc:`datasets/vectors`
-- |-| doc: ``datasets/profile``
 
 v0.1.x Plotting
 ---------------

doc/whatsnew.rst

@@ -29,6 +29,8 @@ v0.2.2 (unreleased)
 New features
 ------------
 
+- Add accessor method :meth:`.Dataset.hyoga.profile`, a new example and a
+  documentation section for profile interpolation (:issue:`18`, :pull:`46`).
 - Add plot method :meth:`.Dataset.hyoga.plot.scale_bar` and a new example for
   automatically sized, anchored scale bars (:issue:`16`, :pull:`44`).
 

examples/interp/plot_profile_altitudes.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+# Copyright (c) 2022, Julien Seguinot (juseg.github.io)
+# GNU General Public License v3.0+ (https://www.gnu.org/licenses/gpl-3.0.txt)
+
+"""
+Profile interpolation
+=====================
+
+Demonstrate interpolating two-dimensional model output onto a one-dimensional
+profile defined by x and y coordinates inside a shapefile.
+"""
+
+import matplotlib.pyplot as plt
+import hyoga
+
+# initialize figure
+fig, (ax, pfax) = plt.subplots(ncols=2, gridspec_kw=dict(width_ratios=(1, 2)))
+
+# open demo data
+with hyoga.open.example('pism.alps.out.2d.nc') as ds:
+
+    # plot 2D model output
+    ds.hyoga.plot.bedrock_altitude(ax=ax, center=False)
+    ds.hyoga.plot.ice_margin(ax=ax, edgecolor='tab:blue', linewidths=1)
+    ds.hyoga.plot.ice_margin(ax=ax, facecolor='tab:blue')
+
+    # interpolate along profile
+    ds = ds.hyoga.profile(hyoga.open.example('profile.rhine.shp'))
+
+    # plot profile line in map view
+    ax.plot(ds.x, ds.y, color='w', dashes=(2, 1))
+    ax.plot(ds.x[0], ds.y[0], 'wo')
+
+    # plot bedrock and surface profiles
+    ds.hyoga.getvar('bedrock_altitude').plot(ax=pfax, color='0.25')
+    ds.hyoga.getvar('surface_altitude').plot(ax=pfax, color='tab:blue')
+
+# set map axes properties
+ax.set_xlim(425e3, 575e3)
+ax.set_ylim(5000e3, 5400e3)
+ax.set_title('map view')
+
+# set profile axes properties
+pfax.set_title('profile view')
+pfax.yaxis.set_label_position("right")
+pfax.yaxis.tick_right()
+
+# show
+plt.show()

hyoga/core/accessor.py

@@ -9,6 +9,7 @@
 """
 
 import warnings
+import geopandas
 import numpy as np
 import scipy.ndimage
 import xarray as xr
@@ -416,6 +417,66 @@ def interp(self, datasource, ax=None, sigma=None, threshold=None):
         # return interpolated data
         return ds
 
+    def profile(self, datasource, interval=None):
+        """Interpolate onto coordinates along a profile.
+
+        Parameters
+        ----------
+        datasource : sequence, array, GeoDataFrame, str, Path or file-like
+            Sequence of (x, y) coordinate tuples, (N, 2) coordinate array,
+            GeoDataFrame or path to a shapefile containing a single line,
+            along which the dataset will be interpolated.
+        interval : float, optional
+            If provided, resample (linearly interpolate) profile coordinates
+            to a fixed spatial resolution given by ``interval``. If ``None``,
+            the data are interpolated to the exact ``datasource`` coordinate,
+            which may produce an irreguar grid.
+
+        Returns
+        -------
+        dataset : Dataset
+            The interpolated dataset, where horizontal dimensions ``x`` and
+            ``y`` are replaced by a new dimension ``d`` with a grid spacing of
+            either ``interval`` or the distance between points in
+            ``datasource``.
+        """
+
+        # read profile from datasource
+        if isinstance(datasource, list):
+            x, y = np.array(datasource).T
+        elif isinstance(datasource, np.ndarray):
+            x, y = datasource.T
+        elif isinstance(datasource, geopandas.GeoDataFrame):
+            x, y = datasource.squeeze().geometry.coords.xy
+        else:
+            datasource = geopandas.read_file(datasource)
+            x, y = datasource.squeeze().geometry.coords.xy
+
+        # compute distance along profile
+        dist = ((np.diff(x)**2+np.diff(y)**2)**0.5).cumsum()
+        dist = np.insert(dist, 0, 0)
+
+        # build coordinate xarrays
+        x = xr.DataArray(x, coords=[dist], dims='d')
+        y = xr.DataArray(y, coords=[dist], dims='d')
+
+        # if interval was given, interpolate coordinates
+        if interval is not None:
+            dist = np.arange(0, dist[-1], interval)
+            x = x.interp(d=dist, method='linear')
+            y = y.interp(d=dist, method='linear')
+
+        # interpolate dataset to new coordinates
+        ds = self._ds.interp(x=x, y=y, method='linear', assume_sorted=True)
+
+        # set new coordinate attributes
+        ds.d.attrs.update(long_name='distance along profile')
+        if 'units' in self._ds.x.attrs:
+            ds.d.attrs.update(units=self._ds.x.units)
+
+        # return interpolated dataset
+        return ds
+
     def where(self, cond, **kwargs):
         """Filter glacier (non-bedrock) variables according to a condition.
 

hyoga/open/example.py

@@ -8,15 +8,32 @@
 series and plotting output from other models.
 """
 
+import os.path
+import geopandas
 import xarray as xr
 
 import hyoga.open.downloader
 
 
 def example(filename='pism.alps.out.2d.nc'):
     """Open cached example dataset from hyoga-data github repository."""
+
+    # github repo url
     repo = 'https://raw.githubusercontent.com/juseg/hyoga-data/main'
-    model = filename.split('.')[0]
-    url = '/'.join((repo, model, filename))
-    path = hyoga.open.downloader.CacheDownloader()(url, filename)
-    return xr.open_dataset(path)
+
+    # open model data with xarray
+    if filename.endswith('.nc'):
+        model = filename.split('.')[0]
+        url = '/'.join((repo, model, filename))
+        path = hyoga.open.downloader.CacheDownloader()(url, filename)
+        return xr.open_dataset(path)
+
+    # open shapefiles with geopandas
+    if filename.endswith('.shp'):
+        url = '/'.join((repo, 'shp', os.path.splitext(filename)[0] + '.zip'))
+        path = os.path.join('examples', 'shp', filename)
+        path = hyoga.open.downloader.ShapeZipDownloader()(url, path, filename)
+        return geopandas.read_file(path)
+
+    # raise error for anything else
+    raise ValueError(f'Could not recognized format for {filename}.')