Fix `BlockwiseCoreg` errors

[rhugonnet]

Since #530, it looks like BlockwiseCoreg has an issue: apply seems wrong, in the documentation built in #502 example it creates more errors than it corrects (wrong direction?). Additionally, it fails with BiasCorr method during the apply step.

[rhugonnet]

Additionally, would be a good moment to add more tests to avoid this happening again in the future.

[rhugonnet]

OK, the cause is indeed that BlockwiseCoreg either applies or estimates the horizontal shift in the wrong direction contrary to classic apply (vertical unaffected), I don't understand why yet...

[rhugonnet]

Even when the coregistration is applied in the right direction (which is not the case by default for BlockwiseCoreg somehow, even after #585, so probably an error in the original implementation), the warp_dem function creates weird artefacts... See example below that I had written in the documentation (temporarily removed for now):

Dividing coregistration in blocks

The {class}~xdem.coreg.BlockwiseCoreg object

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:

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:

# Run 16 block coregistrations
aligned_dem = blockwise.fit_and_apply(ref_dem, tba_dem_shifted)

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

where the tba_dem_shifted is generated above in the coregistration page with:

:tags: [hide-cell]
:mystnb:
:  code_prompt_show: "Show the code for adding a horizontal and vertical shift"
:  code_prompt_hide: "Hide the code for adding a horizontal and vertical shift"

x_shift = 30
y_shift = 30
z_shift = 10
# Affine matrix for 3D transformation
matrix = np.array(
    [
        [1, 0, 0, x_shift],
        [0, 1, 0, y_shift],
        [0, 0, 1, z_shift],
        [0, 0, 0, 1],
    ]
)
# We create misaligned elevation data
tba_dem_shifted = xdem.coreg.apply_matrix(ref_dem, matrix)

And the blockwise before/after plot when I force the direction of shift correction:

We need to find the origin of these artefacts... likely in the warp_dem function.

[erikmannerfelt]

Hmm, do you mean that the artefacts are what can be seen around the edge of the figure that you sent? Meaning that it's not perfect at the edges?

Also, I don't understand how it can suddenly shift the apply in the wrong direction and still produce the right result in the documentation: https://xdem.readthedocs.io/en/stable/advanced_examples/plot_blockwise_coreg.html#sphx-glr-advanced-examples-plot-blockwise-coreg-py

[rhugonnet]

Yes those are the artefacts!

The "stable" documentation is the latest release. The wrong direction compared to other apply appeared since #530, so it's in "latest": https://xdem.readthedocs.io/en/latest/advanced_examples/plot_blockwise_coreg.html