Lack shape consistency on fit and apply

[marinebcht]

Resolves #756

Context

The problem appears when we need to do a fit_and_apply() with a different size mask of the elevation inputs.

  File "....xdem/xdem/coreg/base.py", line 227, in _preprocess_coreg_fit_raster_raster
    invalid_mask = np.logical_or.reduce((~inlier_mask, ref_mask, tba_mask))
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

As a reminder, the fit_and_apply(reference_dem, dem_to_be_aligned) works like this :

ref \ tba	DEM	array	PC
DEM	OK	OK (1/2*) with transform/crs information	OK
array	OK	OK (3*) with transform/crs information	OK (1*) with transform/crs information
PC	OK	OK (1*) with transform/crs information	KO

(1*) error in _preprocess_coreg_apply if not transform/crs

(2*) /!\ warning f"'{name}' of type {type(dem)} overrides the given 'transform'") dans _preprocess_coreg_fit_raster_raster car If any input is a Raster, use its transform if 'transform is None'.

(3*) error in _preprocess_coreg_fit if not transform/crs

Resolution

When the mask is not the same size as elevation data, here are the processes :

mask is a raster

ref \ tba	DEM	array	PC
DEM	reproject with ref (a)	reproject with ref/transform/crs (a)	reproject with ref (b)
array	reproject with tba (a)	reproject with ref/transform/crs (a)	reproject with array/transform/crs (b)
PC	reproject with tba (b)	reproject with array/transform/crs (b)	not possible

(a) _preprocess_coreg_fit_raster_raster
(b) _preprocess_coreg_fit_raster_point

mask is an array

(E2) = No solution found to reproject mask array
=> Raises error "Input mask array can't be a different size array as input elevation."

Tests

The reference results (dem + shift values) are computed with a cropped mask, reprojected with reference_dem.

dem_aligned_crop = nuth_kaab_crop.fit_and_apply(reference_dem, to_be_aligned_dem, inlier_mask_crop_proj, random_state=42)

All nuth_kaab_crop.fit_and_apply with the different configs listed in the previous table + the cropped mask inlier_mask_crop_proj need to match the reference results.

[marinebcht]

I have a problem that I don't understand why it isn't working.

I would like to compare the NuthKaab.fit_and_apply() with :

• large input dems with a cropped mask (N*M)
• all input cropped dems (N*M) and the cropped mask (N*M)

I don't find the same output shifts and I don't know if it's normal or not. I don't think so but I don't understand why my additional code can lead to this error.

import geoutils as gu
import xdem

# get data
reference_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
dem_to_be_aligned = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
inlier_mask = ~glacier_outlines.create_mask(reference_dem)

# crop all data (same size)
nrows, ncols = inlier_mask.shape
inlier_mask_crop = inlier_mask.icrop((0, 0, ncols - 100, nrows - 100))
reference_dem_crop = reference_dem.icrop((0, 0, ncols - 100, nrows - 100))
dem_to_be_aligned_crop = dem_to_be_aligned.icrop((0, 0, ncols - 100, nrows - 100))

list_shift = ["shift_x", "shift_y", "shift_z"]

# run large ref, large tba, cropped mask
nuthkaab_1 = xdem.coreg.NuthKaab()
nuthkaab_1.fit_and_apply(reference_dem, dem_to_be_aligned, inlier_mask=inlier_mask_crop, random_state=42)
shifts_large_large_crop = [nuthkaab_1.meta["outputs"]["affine"][k] for k in list_shift]
print ("shift large ref, large tba, cropped mask:", shifts_large_large_crop)

# run cropped ref, cropped tba, cropped mask
nuthkaab_2 = xdem.coreg.NuthKaab()
nuthkaab_2.fit_and_apply(reference_dem_crop, dem_to_be_aligned_crop, inlier_mask=inlier_mask_crop, random_state=42)
shifts_crop_crop_crop = [nuthkaab_2.meta["outputs"]["affine"][k] for k in list_shift]
print ("shift cropped ref, cropped tba, cropped mask:", shifts_crop_crop_crop)

shift large ref, large tba, cropped mask: [9.17693649962902, 2.534747406844209, -2.069996386286391]
shift cropped ref, cropped tba, cropped mask: [9.909596304477367, 2.6832071805679862, -1.9854913010673556]

[marinebcht]

Another strange problem : I don't find the same output shifts between :

• NuthKaab.fit_and_apply(raster, raster)
• NuthKaab.fit_and_apply(pc, raster)
• NuthKaab.fit_and_apply(raster, pc)

Is that normal ?

import xdem

reference_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
reference_dem_pc = reference_dem.to_pointcloud().ds
reference_dem_pc.rename(columns={'b1': 'z'}, inplace=True)

dem_to_be_aligned = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
dem_to_be_aligned_pc = dem_to_be_aligned.to_pointcloud().ds
dem_to_be_aligned_pc.rename(columns={'b1': 'z'}, inplace=True)

list_shift = ["shift_x", "shift_y", "shift_z"]

# raster/raster
nuthkaab_1 = xdem.coreg.NuthKaab()
nuthkaab_1.fit_and_apply(reference_dem, dem_to_be_aligned, random_state=42)
shifts_raster_raster = [nuthkaab_1.meta["outputs"]["affine"][k] for k in list_shift]

# pc/raster
nuthkaab_2 = xdem.coreg.NuthKaab()
nuthkaab_2.fit_and_apply(reference_dem_pc, dem_to_be_aligned, random_state=42)
shifts_pc_raster = [nuthkaab_2.meta["outputs"]["affine"][k] for k in list_shift]

# raster/pc
nuthkaab_3 = xdem.coreg.NuthKaab()
nuthkaab_3.fit_and_apply(reference_dem, dem_to_be_aligned_pc, random_state=42)
shifts_raster_pc = [nuthkaab_3.meta["outputs"]["affine"][k] for k in list_shift]

raster/raster: [8.547442937129617, 1.2731213536054102, -2.4275404619966707]
pc/raster: [8.707025183613736, 1.1431423299389598, -2.43461892048299]
raster/pc: [8.363922242693267, 1.2221001304179122, -2.4504983449720896]

[marinebcht]

Note : need a validation for (E1) and (E2) @belletva @rhugonnet @adehecq
But it can be of another issue :)

[rhugonnet]

@marinebcht Amazing description and summary, I grasped everything immediately. 😉

And great idea for the extra consistency tests! (crop/uncropped and pointcloud/raster)

Here are possible explanations:

• For uncropped/crop: In the first case (inlier_mask is cropped), the random subsample coordinates are computed using random indices of the N*M array. While, in the second case, they are computed using indices of your cropped array (N - 100 x M - 100). So even with random_state=42, the code is selecting different subsamples.
  If you pin subsample=1 to use all samples in both cases, then there should be no more randomness, and your test should yield the exact same result for both cases! (hopefully)

• For pointcloud/raster input and order: Same as above but with points pre-selected differently du to the various types of the datasets. With subsample=1 you should get the same result hopefully. If not, we have things to fix!

For the PR changes:
(E1): If the mask is a raster, we can run mask.reproject(reference_dem, resampling="nearest") using the _reproject() that works with transform/crs input. But we should raise an error if the mask doesn't intersect the original raster: mask.get_projected_footprint(crs).intersects(dem.get_projected_footprint()), also here using the equivalent with transform/crs input.
(E2): OK for me.

[marinebcht]

@rhugonnet

Good idea but I would raise a warning and not an error. If I understand correctly, if inlier_mask does not interesect ref_dem => inlier_mask is like does not exists (as value = True after reprojecion)

About (E1), you confirm that we can't reproject a raster mask with ...

• ref_dem an array
• tba_dem an array
• transfrom and crs info from ref_dem

Plus, add subsample=1 for the uncropped/crop and pointcloud/raster does not give same shifts :(

[rhugonnet]

Quick answers:

• Inlier mask out of bounds: I think our only option is to raise an error: if the user is passing an inlier_mask that does not intersect the DEMs, then all other values are by definition False (outside the mask), so we can't run the coregistration further.
• Reprojection: No, we can reproject, this is what I meant above. Here is the way: ref_rst = Raster.from_array(data=ref_dem, transform=transform, crs=crs); inlier_mask = inlier_mask.reproject(ref_rst),
• New tests differences: Then another difference is introduced somewhere, we need to dig more... We could start by checking the output of preprocess_fit is indeed the same in your 2 cases. There actually aren't many tests for preprocessing functions right now. Amélie had started some here but we never got to merge them: Dask for coregistration #525

[marinebcht]

1/ Inlier mask out of bounds

I don't know if we talk about the same thing but, to me, if a mask does not intersect the dem it currently make the fit_and_apply works like there is no mask, here a little test to explain it :

reference_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
dem_to_be_aligned = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
inlier_mask = ~glacier_outlines.create_mask(reference_dem)

# crop datas
nrows, ncols = inlier_mask.shape
reference_dem = reference_dem.icrop((ncols - 88, nrows - 88, ncols, nrows))
inlier_mask_crop = inlier_mask.icrop((0, 0, ncols - 100, nrows - 100))

print ("Intersection =", inlier_mask_crop.get_footprint_projected(reference_dem.crs).intersects(reference_dem.get_footprint_projected(reference_dem.crs))[0])

=> Intersection = False

res_nuthkaab_with_mask = xdem.coreg.NuthKaab().fit_and_apply(reference_dem, dem_to_be_aligned, inlier_mask=inlier_mask_crop, random_state=42, subsample=1)
res_nuthkaab_without_mask = xdem.coreg.NuthKaab().fit_and_apply(reference_dem, dem_to_be_aligned, random_state=42, subsample=1)
print ("Res stat Coreg with/without mask =", res_nuthkaab_with_mask.get_stats() == res_nuthkaab_without_mask.get_stats())

=> Res stat Coreg with/without mask = True

Plus, there is no actual test about the intersection (even with size ref_dem = size inlier_mask).

2/ Reprojection

OK ! Got it, thanks !

[rhugonnet]

@marinebcht Yes indeed, the behaviour of 1/ is not correct, it would have to raise an error...

I thought an error/warning was raised by reproject() when Rasters did not intersect, but apparently not. I just checked by adding this line to your example above:

inlier_mask_crop.reproject(reference_dem, resampling="nearest")

Does not fail. It runs, but outputs a boolean raster full of 1s... I would have expected a Raster full of nodata values, but because nodata values don't exist for a boolean raster, they must be converted differently during the operation.

We should probably add a new error/warning in reproject()? (I checked the code and couldn't find anything on this)
Related to: GlacioHack/geoutils#444, GlacioHack/geoutils#495

[marinebcht]

Oh ! Ok so... maybe we can treat this particular case (inlier mask that not intersect input ref dem) in a open/new issues :

• to work on the reprojection error
• add inlier mask intersection test with ref dem (casue actually it is never done, the bug I had only comes from the difference of shape)

@adebardo @belletva

[rhugonnet]

@marinebcht Following the meeting, I tried to understand the differences with cropped inlier_mask/ref/tba VS full inlier_mask/ref/tba you described above!

To circumvent the current error in Raster.reproject() for masks, I used an array of the right size for inlier_mask, which thus doesn't need reprojection, by extending it back to the right size with an array full of False values.
I added subsample=1, and we indeed get the exact sample result:

import geoutils as gu
import xdem
import numpy as np
# get data
reference_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
dem_to_be_aligned = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
inlier_mask = ~glacier_outlines.create_mask(reference_dem)
# crop all data (same size)
nrows, ncols = inlier_mask.shape
inlier_mask_crop = inlier_mask.icrop((0, 0, ncols - 100, nrows - 100))
reference_dem_crop = reference_dem.icrop((0, 0, ncols - 100, nrows - 100))
dem_to_be_aligned_crop = dem_to_be_aligned.icrop((0, 0, ncols - 100, nrows - 100))

# ADDED THIS STEP
inlier_mask_crop_ext = np.zeros(dem_to_be_aligned.shape, dtype=bool)
inlier_mask_crop_ext[0:-100, 0:-100] = inlier_mask_crop.data

list_shift = ["shift_x", "shift_y", "shift_z"]
# run large ref, large tba, cropped mask
nuthkaab_1 = xdem.coreg.NuthKaab(subsample=1)
nuthkaab_1.fit_and_apply(reference_dem, dem_to_be_aligned, inlier_mask=inlier_mask_crop_ext.data, random_state=42)
shifts_large_large_crop = [nuthkaab_1.meta["outputs"]["affine"][k] for k in list_shift]
print ("shift large ref, large tba, cropped mask:", shifts_large_large_crop)
# run cropped ref, cropped tba, cropped mask
nuthkaab_2 = xdem.coreg.NuthKaab(subsample=1)
nuthkaab_2.fit_and_apply(reference_dem_crop, dem_to_be_aligned_crop, inlier_mask=inlier_mask_crop.data, random_state=42)
shifts_crop_crop_crop = [nuthkaab_2.meta["outputs"]["affine"][k] for k in list_shift]
print ("shift cropped ref, cropped tba, cropped mask:", shifts_crop_crop_crop)

shift large ref, large tba, cropped mask: [np.float64(9.891881736564647), np.float64(2.6642006810211187), -1.9917773237056622]
shift cropped ref, cropped tba, cropped mask: [np.float64(9.890835298209975), np.float64(2.6639441268127864), -1.9913401502903412]

Turning off subsample=1, I was surprised that the results were still the same! But actually, that's because the subsampling happens on the indices the flattened array, which happen to be the same because, in the example above, you cropped at the end. If we crop in the beginning also, then results differ without subsample=1 (as expected! 😄).
So it all works as intended! We just need to fix the reprojection with the inlier_mask, and we can add this to the tests/, great idea 😉

Next, I'll look at the pointcloud/raster differences!

[rhugonnet]

Alright, I also figured out the differences for point-raster above!

Some coregistration methods are "point-grid" and thus slightly assymetric because they use the grid to derive X/Y gradients to optimize the coregistration algorithms. This includes NuthKaab, LZD and ICP(method="point-to-plane"). All other methods are fully symmetric.
When converting one dataset to a point cloud, the algorithm cannot use it for a gradient, and thus defaults to the other grid for the gradient, and thus yields a different result.

If we test a symmetric method, such as VerticalShift, or IPC(method="point-to-point"), then the result is exactly the same:

import xdem
reference_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
reference_dem = reference_dem.icrop((0, 0, 100, 100))
reference_dem_pc = reference_dem.to_pointcloud().ds
reference_dem_pc.rename(columns={'b1': 'z'}, inplace=True)
dem_to_be_aligned = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
dem_to_be_aligned = dem_to_be_aligned.icrop((0, 0, 100, 100))
dem_to_be_aligned_pc = dem_to_be_aligned.to_pointcloud().ds
dem_to_be_aligned_pc.rename(columns={'b1': 'z'}, inplace=True)
list_shift = ["shift_x", "shift_y", "shift_z"]
# raster/raster
nuthkaab_1 = xdem.coreg.ICP(subsample=1, method="point-to-point")
nuthkaab_1.fit_and_apply(reference_dem, dem_to_be_aligned, random_state=42)
print([nuthkaab_1.meta["outputs"]["affine"][k] for k in list_shift])
# pc/raster
nuthkaab_2 = xdem.coreg.ICP(subsample=1, method="point-to-point")
nuthkaab_2.fit_and_apply(reference_dem_pc, dem_to_be_aligned, random_state=42)
print([nuthkaab_2.meta["outputs"]["affine"][k] for k in list_shift])
# raster/pc
nuthkaab_3 = xdem.coreg.ICP(subsample=1, method="point-to-point")
nuthkaab_3.fit_and_apply(reference_dem, dem_to_be_aligned_pc, random_state=42)
print([nuthkaab_3.meta["outputs"]["affine"][k] for k in list_shift])

[np.float64(0.04527504684651904), np.float64(-0.02808608440218554), np.float64(-4.176503179789753)]
[np.float64(0.04527504684651904), np.float64(-0.02808608440218554), np.float64(-4.176503179789753)]
[np.float64(0.04527504684651904), np.float64(-0.02808608440218554), np.float64(-4.176503179789753)]

So all seems to be working as expected! 😄
We can also include those tests somewhere @marinebcht, could be in a separate PR!

[marinebcht]

@rhugonnet

Problem 1) Oh yes, I undestand now the consequence of the reprojection error in the case of a mask ^^ Top, I will add the test with the new geoutils branch/version !

Problem 2) Ok, I understand now even it was a little bit strange at first. Maybe we can add it somewhere in the doc ? (The fact that we can have different coreg results the a raster/raster vs raster/point cloud for several method) @belletva @adebardo

[marinebcht]

Added #828 to add some tests for the two points I checked :

• cropped inlier_mask/ref/tba VS full inlier_mask/ref/tba
• pointcloud/raster differences (+ add a little line in the doc for non symetrical function)

[rhugonnet]

This file needs to be restored to match main (wrong merge?).
All these dependencies mirror the ones in setup.cfg, for those using a conda install.

[rhugonnet]

Amazing series of tests, thanks a lot 😉

[rhugonnet]

Perfect! See the small issue with the dev-env file.
We can add the new tests in a separate PR following your issue.