Minor improvements to gallery examples

examples/advanced/plot_slope_methods.py

@@ -30,15 +30,7 @@ def plot_attribute(attribute, cmap, label=None, vlim=None):
     else:
         vlims = {}
 
-    plt.imshow(
-        attribute.squeeze(),
-        cmap=cmap,
-        extent=[dem.bounds.left, dem.bounds.right, dem.bounds.bottom, dem.bounds.top],
-        **vlims,
-    )
-    if label is not None:
-        cbar = plt.colorbar()
-        cbar.set_label(label)
+    attribute.plot(cmap=cmap, cbar_title=label)
 
     plt.xticks([])
     plt.yticks([])
@@ -51,14 +43,14 @@ def plot_attribute(attribute, cmap, label=None, vlim=None):
 # Slope with method of `Horn (1981) <http://dx.doi.org/10.1109/PROC.1981.11918>`_ (GDAL default), based on a refined
 # approximation of the gradient  (page 18, bottom left, and pages 20-21).
 
-slope_horn = xdem.terrain.slope(dem.data, resolution=dem.res)
+slope_horn = xdem.terrain.slope(dem)
 
 plot_attribute(slope_horn, "Reds", "Slope of Horn (1981) (°)")
 
 # %%
 # Slope with method of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.3290120107>`_, Equation 13.
 
-slope_zevenberg = xdem.terrain.slope(dem.data, resolution=dem.res, method="ZevenbergThorne")
+slope_zevenberg = xdem.terrain.slope(dem, method="ZevenbergThorne")
 
 plot_attribute(slope_zevenberg, "Reds", "Slope of Zevenberg and Thorne (1987) (°)")
 
@@ -73,7 +65,7 @@ def plot_attribute(attribute, cmap, label=None, vlim=None):
 # The differences are negative, implying that the method of Horn always provides flatter slopes.
 # Additionally, they seem to occur in places of high curvatures. We verify this by plotting the maximum curvature.
 
-maxc = xdem.terrain.maximum_curvature(dem.data, resolution=dem.res)
+maxc = xdem.terrain.maximum_curvature(dem)
 
 plot_attribute(maxc, "RdYlBu", "Maximum curvature (100 m $^{-1}$)", vlim=2)
 
@@ -102,11 +94,11 @@ def plot_attribute(attribute, cmap, label=None, vlim=None):
 # We perform the same exercise to analyze the differences in terrain aspect. We compute the difference modulo 360°,
 # to account for the circularity of aspect.
 
-aspect_horn = xdem.terrain.aspect(dem.data)
-aspect_zevenberg = xdem.terrain.aspect(dem.data, method="ZevenbergThorne")
+aspect_horn = xdem.terrain.aspect(dem)
+aspect_zevenberg = xdem.terrain.aspect(dem, method="ZevenbergThorne")
 
 diff_aspect = aspect_horn - aspect_zevenberg
-diff_aspect_mod = np.minimum(np.mod(diff_aspect, 360), 360 - np.mod(diff_aspect, 360))
+diff_aspect_mod = np.minimum(diff_aspect % 360, 360 - diff_aspect % 360)
 
 plot_attribute(
     diff_aspect_mod, "Spectral", "Aspect of Horn (1981) minus\n aspect of Zevenberg and Thorne (1987) (°)", vlim=[0, 90]

examples/advanced/plot_variogram_estimation_modelling.py

@@ -78,7 +78,9 @@
 # conveniently by :func:`xdem.spatialstats.sample_empirical_variogram`:
 # Dowd's variogram is used for robustness in conjunction with the NMAD (see :ref:`robuststats-corr`).
 
-df = xdem.spatialstats.sample_empirical_variogram(values=dh, subsample=1000, n_variograms=10, estimator="dowd", random_state=42)
+df = xdem.spatialstats.sample_empirical_variogram(
+    values=dh, subsample=1000, n_variograms=10, estimator="dowd", random_state=42
+)
 
 # %%
 # *Note: in this example, we add a* ``random_state`` *argument to yield a reproducible random sampling of pixels within

examples/basic/plot_dem_subtraction.py

@@ -1,80 +1,60 @@
 """
-DEM subtraction
-===============
+DEM differencing
+================
 
-Subtracting one DEM with another should be easy!
-This is why ``xdem`` (with functionality from `geoutils <https://geoutils.readthedocs.io>`_) allows directly using the ``-`` or ``+`` operators on :class:`xdem.DEM` objects.
+Subtracting a DEM with another one should be easy.
 
-Before DEMs can be compared, they need to be reprojected/resampled/cropped to fit the same grid.
-The :func:`xdem.DEM.reproject` method takes care of this.
+xDEM allows to use any operator on :class:`xdem.DEM` objects, such as :func:`+<operator.add>` or :func:`-<operator.sub>` as well as most NumPy functions
+while respecting nodata values and checking that georeferencing is consistent. This functionality is inherited from `GeoUtils' Raster class <https://geoutils.readthedocs.io>`_.
+
+Before DEMs can be compared, they need to be reprojected to the same grid and have the same 3D CRSs. The :func:`geoutils.Raster.reproject` and :func:`xdem.DEM.to_vcrs` methods are used for this.
 
 """
 import geoutils as gu
-import matplotlib.pyplot as plt
 
 import xdem
 
 # %%
+# We load two DEMs near Longyearbyen, Svalbard.
 
 dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
 dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem_coreg"))
 
 # %%
-# We can print the information about the DEMs for a "sanity check"
+# We can print the information about the DEMs for a "sanity check".
 
-print(dem_2009)
-print(dem_1990)
+dem_2009.info()
+dem_1990.info()
 
 # %%
 # In this particular case, the two DEMs are already on the same grid (they have the same bounds, resolution and coordinate system).
-# If they don't, we need to reproject one DEM to fit the other.
-# :func:`xdem.DEM.reproject` is a multi-purpose method that ensures a fit each time:
+# If they don't, we need to reproject one DEM to fit the other using :func:`geoutils.Raster.reproject`:
 
-_ = dem_1990.reproject(dem_2009)
+dem_1990 = dem_1990.reproject(dem_2009)
 
 # %%
 # Oops!
-# ``xdem`` just warned us that ``dem_1990`` did not need reprojection, but we asked it to anyway.
-# To hide this prompt, add ``.reproject(..., silent=True)``.
+# GeoUtils just warned us that ``dem_1990`` did not need reprojection. We can hide this output with ``.reproject(..., silent=True)``.
 # By default, :func:`xdem.DEM.reproject` uses "bilinear" resampling (assuming resampling is needed).
-# Other options are "nearest" (fast but inaccurate), "cubic_spline", "lanczos" and others.
-# See `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.raster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_ for more information about reprojection.
+# Other options are detailed at `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.raster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_.
 #
-# Now, we are ready to generate the dDEM:
+# We now compute the difference by simply substracting, passing `stats=True` to :func:`geoutils.Raster.info` to print statistics.
 
 ddem = dem_2009 - dem_1990
 
-print(ddem)
+ddem.info(stats=True)
 
 # %%
 # It is a new :class:`xdem.DEM` instance, loaded in memory.
-# Let's visualize it:
-
-ddem.plot(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
-
-# %%
-# Let's add some glacier outlines
+# Let's visualize it, with some glacier outlines.
 
 # Load the outlines
 glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
-
-# Need to create a common matplotlib Axes to plot both on the same figure
-# The xlim/ylim commands are necessary only because outlines extend further than the raster extent
-ax = plt.subplot(111)
-ddem.plot(ax=ax, cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
-glacier_outlines.ds.plot(ax=ax, fc="none", ec="k")
-plt.xlim(ddem.bounds.left, ddem.bounds.right)
-plt.ylim(ddem.bounds.bottom, ddem.bounds.top)
-plt.title("With glacier outlines")
-plt.show()
-
-# %%
-# For missing values, ``xdem`` provides a number of interpolation methods which are shown in the other examples.
+glacier_outlines = glacier_outlines.crop(ddem, clip=True)
+ddem.plot(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
+glacier_outlines.plot(ref_crs=ddem, fc="none", ec="k")
 
 # %%
-# Saving the output to a file is also very simple
+# And we save the output to file.
 
 ddem.save("temp.tif")
-
-# %%
-# ... and that's it!

examples/basic/plot_icp_coregistration.py

@@ -1,12 +1,12 @@
 """
-Iterative Closest Point coregistration
+Iterative closest point coregistration
 ======================================
-Some DEMs may for one or more reason be erroneously rotated in the X, Y or Z directions.
-Established coregistration approaches like :ref:`coregistration-nuthkaab` work great for X, Y and Z *translations*, but rotation is not accounted for at all.
 
-Iterative Closest Point (ICP) is one method that takes both rotation and translation into account.
-It is however not as good as :ref:`coregistration-nuthkaab` when it comes to sub-pixel accuracy.
-Fortunately, ``xdem`` provides the best of two worlds by allowing a combination of the two.
+Iterative Closest Point (ICP) is a registration methods accounting for both rotation and translation.
+
+It is used primarily to correct rotations, as it performs worse than :ref:`coregistration-nuthkaab` for sub-pixel shifts.
+Fortunately, xDEM provides the best of two worlds by allowing a combination of the two methods in a pipeline,
+demonstrated below!
 
 **Reference**: `Besl and McKay (1992) <https://doi.org/10.1117/12.57955>`_.
 """
@@ -17,7 +17,7 @@
 import xdem
 
 # %%
-# Let's load a DEM and crop it to a single mountain on Svalbard, called Battfjellet.
+# We load a DEM and crop it to a single mountain on Svalbard, called Battfjellet.
 # Its aspects vary in every direction, and is therefore a good candidate for coregistration exercises.
 dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))