Skip to content

Commit

Permalink
Minor improvements to gallery examples (GlacioHack#493)
Browse files Browse the repository at this point in the history
  • Loading branch information
rhugonnet authored Mar 20, 2024
1 parent 1ad4ad1 commit 48bbf17
Show file tree
Hide file tree
Showing 4 changed files with 39 additions and 65 deletions.
22 changes: 7 additions & 15 deletions examples/advanced/plot_slope_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -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([])
Expand All @@ -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) (°)")

Expand All @@ -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)

Expand Down Expand Up @@ -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]
Expand Down
4 changes: 3 additions & 1 deletion examples/advanced/plot_variogram_estimation_modelling.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
64 changes: 22 additions & 42 deletions examples/basic/plot_dem_subtraction.py
Original file line number Diff line number Diff line change
@@ -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!
14 changes: 7 additions & 7 deletions examples/basic/plot_icp_coregistration.py
Original file line number Diff line number Diff line change
@@ -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>`_.
"""
Expand All @@ -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"))

Expand Down

0 comments on commit 48bbf17

Please sign in to comment.