diff --git a/examples/advanced/plot_slope_methods.py b/examples/advanced/plot_slope_methods.py index 9eb4063e..5c005c97 100644 --- a/examples/advanced/plot_slope_methods.py +++ b/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) `_ (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) `_, 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] diff --git a/examples/advanced/plot_variogram_estimation_modelling.py b/examples/advanced/plot_variogram_estimation_modelling.py index 3a03deea..43cff5a6 100644 --- a/examples/advanced/plot_variogram_estimation_modelling.py +++ b/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 diff --git a/examples/basic/plot_dem_subtraction.py b/examples/basic/plot_dem_subtraction.py index f5def25f..ee158338 100644 --- a/examples/basic/plot_dem_subtraction.py +++ b/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 `_) 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:`+` or :func:`-` as well as most NumPy functions +while respecting nodata values and checking that georeferencing is consistent. This functionality is inherited from `GeoUtils' Raster class `_. + +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() `_ and `rasterio.enums.Resampling `_ for more information about reprojection. +# Other options are detailed at `geoutils.Raster.reproject() `_ and `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! diff --git a/examples/basic/plot_icp_coregistration.py b/examples/basic/plot_icp_coregistration.py index 2fbf28e6..f26eef70 100644 --- a/examples/basic/plot_icp_coregistration.py +++ b/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) `_. """ @@ -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"))