Why can coregistration result in the appearance of block-like (or patchy) artifacts?

[TreeYu123]

I did not use xdem.coreg.BlockwiseCoreg, but instead applied xdem.coreg.NuthKaab; this blocky effect was not present in the original DEMs.

[adebardo]

Thank you for your question. We'll do our best to answer it, could we get a bit more context?
Such as the input DEMs? Or perhaps your script?

Otherwise, off the top of my head, I suspect a very large difference in resolution.

Thanks in advance.

[rhugonnet]

Weird artefacts indeed!

To give more context on how this works under-the-hood in xDEM: Once the shifts are estimated by Coreg.fit(), the apply_matrix function (https://xdem.readthedocs.io/en/stable/gen_modules/xdem.coreg.apply_matrix.html#xdem.coreg.apply_matrix) is applied to the DEM during Coreg.apply(). This apply step is the only one modifying values here.

In the case of Nuth and Kääb, because the matrix is a simple shift (no rotations or scaling), the apply_matrix step is equivalent to a reprojection in the same CRS with reproject().

The apply matrix for all steps is tested here:

xdem/tests/test_coreg/test_base.py

Line 897 in ab155e9

def test_apply_matrix__raster(self, regrid_method: None | str, matrix: NDArrayf) -> None:

And the shift part of it here:

xdem/tests/test_coreg/test_affine.py

Line 84 in ab155e9

def test_reproject_horizontal_shift_samecrs__gdal(self, xoff_yoff: tuple[float, float], get_test_data_path) -> None:

These tests show that outputs agree exactly or within a very small tolerance with other transformations (Shapely at point scale, and GDAL reprojection at grid-scale), so it would be quite surprising if this was an issue from apply_matrix...

But I think I have seen this type of artefact before, long ago. Might have been from the GDAL reprojection. As @adebardo says, we need more context on your code/data to know!

(P.S: Moved this issue to discussions)

[TreeYu123]

Thank you so much. I provided details of two DEMs and codes for your information.

The projection info. of KH-9 DEM

The projection info. of SRTM DEM

Codes:

ref_dem = xdem.DEM(srtm00_dem)
align_kh = xdem.DEM('kh9.tif').reproject(ref_dem)
gla=gu.Vector(pth+'\\MGI\\MGI_1974.shp')
inlier_mask = ~gla.create_mask(ref_dem)
nuth_kaab = xdem.coreg.NuthKaab()
aligned_kh = nuth_kaab.fit_and_apply(ref_dem, align_aster, inlier_mask)
dh79_00 = ref_dem - aligned_kh
dh79_00.save(pthout+'\\kh_srtm_dh.tif')

[rhugonnet]

Hi @TreeYu123,
Sorry we can't work off this. We need either the files of the DEMs (SRTM and KH9 are both open?), or synthetic data reproducing the same artefacts in a minimal reproducible example.