|
1 | 1 | """ |
2 | | -DEM subtraction |
3 | | -=============== |
| 2 | +DEM differencing |
| 3 | +================ |
4 | 4 |
|
5 | | -Subtracting one DEM with another should be easy! |
6 | | -This is why ``xdem`` (with functionality from `geoutils <https://geoutils.readthedocs.io>`_) allows directly using the ``-`` or ``+`` operators on :class:`xdem.DEM` objects. |
| 5 | +Subtracting a DEM with another one should be easy. |
7 | 6 |
|
8 | | -Before DEMs can be compared, they need to be reprojected/resampled/cropped to fit the same grid. |
9 | | -The :func:`xdem.DEM.reproject` method takes care of this. |
| 7 | +xDEM allows directly using any operator on :class:`xdem.DEM` objects, such as :func:`+<operator.add>` or :func:``-<operator.sub>`` as well as most NumPy functions, |
| 8 | +while respecting nodata values and checking that georeferencing is consistent. |
| 9 | +This functionality is inherited from `GeoUtils' Raster class <https://geoutils.readthedocs.io>`_. |
| 10 | +
|
| 11 | +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. |
10 | 12 |
|
11 | 13 | """ |
12 | 14 | import geoutils as gu |
13 | | -import matplotlib.pyplot as plt |
14 | | - |
15 | 15 | import xdem |
16 | 16 |
|
17 | 17 | # %% |
| 18 | +# We load two DEMs near Longyearbyen, Svalbard. |
18 | 19 |
|
19 | 20 | dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem")) |
20 | 21 | dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem_coreg")) |
21 | 22 |
|
22 | 23 | # %% |
23 | | -# We can print the information about the DEMs for a "sanity check" |
| 24 | +# We can print the information about the DEMs for a "sanity check". |
24 | 25 |
|
25 | | -print(dem_2009) |
26 | | -print(dem_1990) |
| 26 | +dem_2009.info() |
| 27 | +dem_1990.info() |
27 | 28 |
|
28 | 29 | # %% |
29 | 30 | # In this particular case, the two DEMs are already on the same grid (they have the same bounds, resolution and coordinate system). |
30 | | -# If they don't, we need to reproject one DEM to fit the other. |
31 | | -# :func:`xdem.DEM.reproject` is a multi-purpose method that ensures a fit each time: |
| 31 | +# If they don't, we need to reproject one DEM to fit the other using :func:`geoutils.Raster.reproject`: |
32 | 32 |
|
33 | | -_ = dem_1990.reproject(dem_2009) |
| 33 | +dem_1990 = dem_1990.reproject(dem_2009) |
34 | 34 |
|
35 | 35 | # %% |
36 | 36 | # Oops! |
37 | | -# ``xdem`` just warned us that ``dem_1990`` did not need reprojection, but we asked it to anyway. |
38 | | -# To hide this prompt, add ``.reproject(..., silent=True)``. |
| 37 | +# GeoUtils just warned us that ``dem_1990`` did not need reprojection. We can hide this output with ``.reproject(..., silent=True)``. |
39 | 38 | # By default, :func:`xdem.DEM.reproject` uses "bilinear" resampling (assuming resampling is needed). |
40 | | -# Other options are "nearest" (fast but inaccurate), "cubic_spline", "lanczos" and others. |
41 | | -# 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. |
| 39 | +# 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>`_. |
42 | 40 | # |
43 | | -# Now, we are ready to generate the dDEM: |
| 41 | +# We now compute the difference by simply substracting, passing `stats=True` to :func:`geoutils.Raster.info` to print statistics. |
44 | 42 |
|
45 | 43 | ddem = dem_2009 - dem_1990 |
46 | 44 |
|
47 | | -print(ddem) |
| 45 | +ddem.info(stats=True) |
48 | 46 |
|
49 | 47 | # %% |
50 | 48 | # It is a new :class:`xdem.DEM` instance, loaded in memory. |
51 | | -# Let's visualize it: |
52 | | - |
53 | | -ddem.plot(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)") |
54 | | - |
55 | | -# %% |
56 | | -# Let's add some glacier outlines |
| 49 | +# Let's visualize it, with some glacier outlines. |
57 | 50 |
|
58 | 51 | # Load the outlines |
59 | 52 | glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines")) |
60 | | - |
61 | | -# Need to create a common matplotlib Axes to plot both on the same figure |
62 | | -# The xlim/ylim commands are necessary only because outlines extend further than the raster extent |
63 | | -ax = plt.subplot(111) |
64 | | -ddem.plot(ax=ax, cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)") |
65 | | -glacier_outlines.ds.plot(ax=ax, fc="none", ec="k") |
66 | | -plt.xlim(ddem.bounds.left, ddem.bounds.right) |
67 | | -plt.ylim(ddem.bounds.bottom, ddem.bounds.top) |
68 | | -plt.title("With glacier outlines") |
69 | | -plt.show() |
70 | | - |
71 | | -# %% |
72 | | -# For missing values, ``xdem`` provides a number of interpolation methods which are shown in the other examples. |
| 53 | +ddem.plot(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)") |
| 54 | +glacier_outlines.plot(ref_crs=ddem, fc="none", ec="k") |
73 | 55 |
|
74 | 56 | # %% |
75 | | -# Saving the output to a file is also very simple |
| 57 | +# And we save the output to file. |
76 | 58 |
|
77 | 59 | ddem.save("temp.tif") |
78 | | - |
79 | | -# %% |
80 | | -# ... and that's it! |
0 commit comments