Fix BlockwiseCoreg shift direction (#665)

.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/

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),

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([])
+```

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)")

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

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,