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

[erikmannerfelt]

Oof yeah that's bad..

It would be interesting to know what the blockwise.stats() would look like before/after the update. That would narrow the problem down a bit.

I've looked through the code and can't find an obvious place where this issue arises. Any ideas?

[rhugonnet]

I could not find much either. Here's what I put together:

• For the wrong direction: it seems directly linked to the way apply is done in BlockwiseCoreg, as all other Coreg functions are now robustly tested for the direction of correction (since Homogenize sign of translation coregistration and add affine tests #585) and use the same apply_matrix. The error probably didn't appear before because NuthKaab had its own _apply_rst function which was overridding the apply_matrix use (and was likely inconsistent with the normal direction of apply). Did we ever test BlockwiseCoreg with anything else than NuthKaab?
• For the artefacts: it seems linked to the way apply + warp_dem work, which I have trouble understanding in detail. If that's OK, I let you dig into this as you know the code better 😉.

[erikmannerfelt]

@rhugonnet, sorry I'm doing my best to keep up! Were these direction issues fixed with the latest PRs or is it still a problem?

Regarding the artefacts, I suspect that it's not warp_dem that's at fault but how the co-registration is extrapolated at the edges. Not sure but that's my suspicion. I can try to make time to figure it out, but I would consider this to be a separate issue altogether!

[rhugonnet]

Still not fixed but will maybe be looked at in #598? @adebardo can keep us posted on this (and potentially warp_dem too).