From 8498fc42d25dfd9d52cd6f536c80c8b1e0f370d7 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Thu, 5 Dec 2024 18:11:30 -0900 Subject: [PATCH] Fix `BlockwiseCoreg` shift direction (#665) --- .gitignore | 3 ++- doc/source/conf.py | 3 +-- doc/source/coregistration.md | 41 +++++++++++++++++++++++++++++++++++ doc/source/uncertainty.md | 2 +- tests/test_coreg/test_base.py | 2 +- xdem/coreg/base.py | 6 ++--- 6 files changed, 49 insertions(+), 8 deletions(-) diff --git a/.gitignore b/.gitignore index 9d8091fb..10b3931c 100644 --- a/.gitignore +++ b/.gitignore @@ -153,5 +153,6 @@ doc/source/gen_modules/ doc/source/sg_execution_times.rst examples/basic/temp.tif -# Directory where myst_nb executes jupyter code +# Directory where myst_nb executes jupyter code and cache doc/jupyter_execute/ +doc/.jupyter_cache/ diff --git a/doc/source/conf.py b/doc/source/conf.py index c8d5dd7a..f37be141 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -66,8 +66,7 @@ nb_execution_raise_on_error = True # To fail documentation build on notebook execution error nb_execution_show_tb = True # To show full traceback on notebook execution error nb_output_stderr = "warn" # To warn if an error is raised in a notebook cell (if intended, override to "show" in cell) - -# autosummary_generate = True +nb_execution_mode = "cache" intersphinx_mapping = { "python": ("https://docs.python.org/", None), diff --git a/doc/source/coregistration.md b/doc/source/coregistration.md index d0b6119c..218fe7f6 100644 --- a/doc/source/coregistration.md +++ b/doc/source/coregistration.md @@ -493,3 +493,44 @@ These metadata are only inputs specific to a given method, outlined in the metho For instance, for {class}`xdem.coreg.Deramp`, an input `poly_order` to define the polynomial order used for the fit, and for {class}`xdem.coreg.DirectionalBias`, an input `angle` to define the angle at which to do the directional correction. + +## Dividing coregistration in blocks + +### The {class}`~xdem.coreg.BlockwiseCoreg` object + +```{caution} +The {class}`~xdem.coreg.BlockwiseCoreg` feature is still experimental: it might not support all coregistration +methods, and create edge artefacts. +``` + +Sometimes, we want to split a coregistration across different spatial subsets of an elevation dataset, running that +method independently in each subset. A {class}`~xdem.coreg.BlockwiseCoreg` can be constructed for this: + +```{code-cell} ipython3 +blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), subdivision=16) +``` + +The subdivision corresponds to an equal-length block division across the extent of the elevation dataset. It needs +to be a number of the form 2{sup}`n` (such as 4 or 256). + +It is run the same way as other coregistrations: + +```{code-cell} ipython3 +# Run 16 block coregistrations +aligned_dem = blockwise.fit_and_apply(ref_dem, tba_dem_shifted) +``` + +```{code-cell} ipython3 +:tags: [hide-input] +:mystnb: +: code_prompt_show: "Show plotting code" +: code_prompt_hide: "Hide plotting code" + +# Plot before and after +f, ax = plt.subplots(1, 2) +ax[0].set_title("Before block NK") +(tba_dem_shifted - ref_dem).plot(cmap='RdYlBu', vmin=-30, vmax=30, ax=ax[0]) +ax[1].set_title("After block NK") +(aligned_dem - ref_dem).plot(cmap='RdYlBu', vmin=-30, vmax=30, ax=ax[1], cbar_title="Elevation differences (m)") +_ = ax[1].set_yticklabels([]) +``` diff --git a/doc/source/uncertainty.md b/doc/source/uncertainty.md index 095149b6..f30e0855 100644 --- a/doc/source/uncertainty.md +++ b/doc/source/uncertainty.md @@ -79,7 +79,7 @@ tba_dem_coreg = tba_dem.coregister_3d(ref_dem, inlier_mask=inlier_mask) ```{code-cell} ipython3 # Estimate elevation uncertainty assuming both DEMs have similar precision -sig_dem, rho_sig = tba_dem_coreg.estimate_uncertainty(ref_dem, stable_terrain=inlier_mask, precision_of_other="same") +sig_dem, rho_sig = tba_dem_coreg.estimate_uncertainty(ref_dem, stable_terrain=inlier_mask, precision_of_other="same", random_state=42) # The error map variability is estimated from slope and curvature by default sig_dem.plot(cmap="Purples", cbar_title=r"Error in elevation (1$\sigma$, m)") diff --git a/tests/test_coreg/test_base.py b/tests/test_coreg/test_base.py index 34f8529f..a6d3d976 100644 --- a/tests/test_coreg/test_base.py +++ b/tests/test_coreg/test_base.py @@ -1150,7 +1150,7 @@ def test_warp_dem() -> None: ) # The warped DEM should have the value 'elev_shift' in the upper left corner. - assert warped_dem[0, 0] == elev_shift + assert warped_dem[0, 0] == -elev_shift # The corner should be zero, so the corner pixel (represents the corner minus resolution / 2) should be close. # We select the pixel before the corner (-2 in X-axis) to avoid the NaN propagation on the bottom row. assert warped_dem[-2, -1] < 1 diff --git a/xdem/coreg/base.py b/xdem/coreg/base.py index 242278e3..5222ddc3 100644 --- a/xdem/coreg/base.py +++ b/xdem/coreg/base.py @@ -3393,8 +3393,8 @@ def _apply_rst( warped_dem = warp_dem( dem=elev, transform=transform, - source_coords=all_points[:, :, 0], - destination_coords=all_points[:, :, 1], + source_coords=all_points[:, :, 1], + destination_coords=all_points[:, :, 0], resampling="linear", ) @@ -3538,7 +3538,7 @@ def warp_dem( if not no_vertical: grid_offsets = scipy.interpolate.griddata( points=destination_coords_scaled[:, :2], - values=destination_coords_scaled[:, 2] - source_coords_scaled[:, 2], + values=source_coords_scaled[:, 2] - destination_coords_scaled[:, 2], xi=(grid_x, grid_y), method=resampling, fill_value=np.nan,