Skip to content

Commit

Permalink
Make coregistration metadata consistent and public (#526)
Browse files Browse the repository at this point in the history
  • Loading branch information
rhugonnet committed May 23, 2024
1 parent 2c9f2b3 commit c30dbec
Show file tree
Hide file tree
Showing 8 changed files with 219 additions and 195 deletions.
2 changes: 1 addition & 1 deletion doc/source/coregistration.md
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,7 @@ vshift.fit(ref_dem, tba_dem, inlier_mask=inlier_mask)
shifted_dem = vshift.apply(tba_dem)
# Use median shift instead
vshift_median = coreg.VerticalShift(vshift_func=np.median)
vshift_median = coreg.VerticalShift(vshift_reduc_func=np.median)
```

## ICP
Expand Down
53 changes: 29 additions & 24 deletions tests/test_coreg/test_affine.py
Original file line number Diff line number Diff line change
Expand Up @@ -74,19 +74,21 @@ def test_vertical_shift(self) -> None:
# Fit the vertical shift model to the data
vshiftcorr.fit(**self.fit_params)

res = self.ref.res[0]

# Check that a vertical shift was found.
assert vshiftcorr._meta.get("vshift") is not None
assert vshiftcorr._meta["vshift"] != 0.0
assert vshiftcorr.meta.get("shift_z") is not None
assert vshiftcorr.meta["shift_z"] != 0.0

# Copy the vertical shift to see if it changes in the test (it shouldn't)
vshift = copy.copy(vshiftcorr._meta["vshift"])
vshift = copy.copy(vshiftcorr.meta["shift_z"])

# Check that the to_matrix function works as it should
matrix = vshiftcorr.to_matrix()
assert matrix[2, 3] == vshift, matrix

# Check that the first z coordinate is now the vertical shift
assert all(vshiftcorr.apply(self.points)["z"].values == vshiftcorr._meta["vshift"])
assert all(vshiftcorr.apply(self.points)["z"].values == vshiftcorr.meta["shift_z"])

# Apply the model to correct the DEM
tba_unshifted, _ = vshiftcorr.apply(self.tba.data, transform=self.ref.transform, crs=self.ref.crs)
Expand All @@ -104,13 +106,13 @@ def test_vertical_shift(self) -> None:
inlier_mask=self.inlier_mask,
)
# Test the vertical shift
newmeta: CoregDict = vshiftcorr2._meta
new_vshift = newmeta["vshift"]
assert np.abs(new_vshift) < 0.01
newmeta: CoregDict = vshiftcorr2.meta
new_vshift = newmeta["shift_z"]
assert np.abs(new_vshift) * res < 0.01

# Check that the original model's vertical shift has not changed
# (that the _meta dicts are two different objects)
assert vshiftcorr._meta["vshift"] == vshift
# (that the _.meta dicts are two different objects)
assert vshiftcorr.meta["shift_z"] == vshift

def test_all_nans(self) -> None:
"""Check that the coregistration approaches fail gracefully when given only nans."""
Expand Down Expand Up @@ -141,9 +143,10 @@ def test_coreg_example(self, verbose: bool = False) -> None:
nuth_kaab = xdem.coreg.NuthKaab()
nuth_kaab.fit(self.ref, self.tba, inlier_mask=self.inlier_mask, verbose=verbose, random_state=42)

# Check the output metadata is always the same
shifts = (nuth_kaab._meta["offset_east_px"], nuth_kaab._meta["offset_north_px"], nuth_kaab._meta["vshift"])
assert shifts == pytest.approx((-0.463, -0.1339999, -1.9922009))
# Check the output .metadata is always the same
shifts = (nuth_kaab.meta["shift_x"], nuth_kaab.meta["shift_y"], nuth_kaab.meta["shift_z"])
res = self.ref.res[0]
assert shifts == pytest.approx((-0.463 * res, -0.1339999 * res, -1.9922009))

def test_gradientdescending(self, subsample: int = 10000, inlier_mask: bool = True, verbose: bool = False) -> None:
"""
Expand All @@ -165,8 +168,9 @@ def test_gradientdescending(self, subsample: int = 10000, inlier_mask: bool = Tr
random_state=42,
)

shifts = (gds._meta["offset_east_px"], gds._meta["offset_north_px"], gds._meta["vshift"])
assert shifts == pytest.approx((0.03525, -0.59775, -2.39144), abs=10e-5)
res = self.ref.res[0]
shifts = (gds.meta["shift_x"], gds.meta["shift_y"], gds.meta["shift_z"])
assert shifts == pytest.approx((0.03525 * res, -0.59775 * res, -2.39144), abs=10e-5)

@pytest.mark.parametrize("shift_px", [(1, 1), (2, 2)]) # type: ignore
@pytest.mark.parametrize("coreg_class", [coreg.NuthKaab, coreg.GradientDescending, coreg.ICP]) # type: ignore
Expand Down Expand Up @@ -205,19 +209,19 @@ def test_coreg_example_shift(self, shift_px, coreg_class, points_or_raster, verb
# The ICP fit only creates a matrix and doesn't normally show the alignment in pixels
# Since the test is formed to validate pixel shifts, these calls extract the approximate pixel shift
# from the matrix (it's not perfect since rotation/scale can change it).
coreg_obj._meta["offset_east_px"] = -matrix[0][3] / res
coreg_obj._meta["offset_north_px"] = -matrix[1][3] / res
coreg_obj.meta["shift_x"] = -matrix[0][3]
coreg_obj.meta["shift_y"] = -matrix[1][3]

# ICP can never be expected to be much better than 1px on structured data, as its implementation often finds a
# minimum between two grid points. This is clearly warned for in the documentation.
precision = 1e-2 if coreg_class.__name__ != "ICP" else 1

if coreg_obj._meta["offset_east_px"] == pytest.approx(-shift_px[0], rel=precision) and coreg_obj._meta[
"offset_north_px"
] == pytest.approx(-shift_px[0], rel=precision):
if coreg_obj.meta["shift_x"] == pytest.approx(-shift_px[0] * res, rel=precision) and coreg_obj.meta[
"shift_y"
] == pytest.approx(-shift_px[0] * res, rel=precision):
return
best_east_diff = coreg_obj._meta["offset_east_px"] - shift_px[0]
best_north_diff = coreg_obj._meta["offset_north_px"] - shift_px[1]
best_east_diff = coreg_obj.meta["shift_x"] - shift_px[0]
best_north_diff = coreg_obj.meta["shift_y"] - shift_px[1]

raise AssertionError(f"Diffs are too big. east: {best_east_diff:.2f} px, north: {best_north_diff:.2f} px")

Expand All @@ -243,9 +247,10 @@ def test_nuth_kaab(self) -> None:
)

# Make sure that the estimated offsets are similar to what was synthesized.
assert nuth_kaab._meta["offset_east_px"] == pytest.approx(pixel_shift, abs=0.03)
assert nuth_kaab._meta["offset_north_px"] == pytest.approx(0, abs=0.03)
assert nuth_kaab._meta["vshift"] == pytest.approx(-vshift, 0.03)
res = self.ref.res[0]
assert nuth_kaab.meta["shift_x"] == pytest.approx(pixel_shift * res, abs=0.03)
assert nuth_kaab.meta["shift_y"] == pytest.approx(0, abs=0.03)
assert nuth_kaab.meta["shift_z"] == pytest.approx(-vshift, 0.03)

# Apply the estimated shift to "revert the DEM" to its original state.
unshifted_dem, _ = nuth_kaab.apply(shifted_dem, transform=self.ref.transform, crs=self.ref.crs)
Expand Down
79 changes: 38 additions & 41 deletions tests/test_coreg/test_base.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,12 +67,11 @@ def test_copy(self, coreg_class: Callable[[], Coreg]) -> None:
corr = coreg_class()
corr_copy = corr.copy()

# Assign some attributes and metadata after copying, respecting the CoregDict type class
corr.vshift = 1
corr._meta["resolution"] = 30
# Assign some attributes and .metadata after copying, respecting the CoregDict type class
corr._meta["shift_z"] = 30
# Make sure these don't appear in the copy
assert corr_copy._meta != corr._meta
assert not hasattr(corr_copy, "vshift")
assert corr_copy.meta != corr.meta
assert not hasattr(corr_copy, "shift_z")

def test_error_method(self) -> None:
"""Test different error measures."""
Expand All @@ -90,7 +89,7 @@ def test_error_method(self) -> None:
assert vshiftcorr.error(dem1, dem2, transform=affine, crs=crs, error_type="median") == 0

# Remove the vertical shift fit and see what happens.
vshiftcorr._meta["vshift"] = 0
vshiftcorr.meta["shift_z"] = 0
# Now it should be equal to dem1 - dem2
assert vshiftcorr.error(dem1, dem2, transform=affine, crs=crs, error_type="median") == -2

Expand All @@ -108,7 +107,7 @@ def test_get_subsample_on_valid_mask(self, subsample: float | int) -> None:
rng = np.random.default_rng(42)
valid_mask = rng.integers(low=0, high=2, size=(width, height), dtype=bool)

# Define a class with a subsample and random_state in the metadata
# Define a class with a subsample and random_state in the .metadata
coreg = Coreg(meta={"subsample": subsample, "random_state": 42})
subsample_mask = coreg._get_subsample_on_valid_mask(valid_mask=valid_mask)

Expand Down Expand Up @@ -141,32 +140,32 @@ def test_subsample(self, coreg: Callable) -> None: # type: ignore
# Check that default value is set properly
coreg_full = coreg()
argspec = inspect.getfullargspec(coreg)
assert coreg_full._meta["subsample"] == argspec.defaults[argspec.args.index("subsample") - 1] # type: ignore
assert coreg_full.meta["subsample"] == argspec.defaults[argspec.args.index("subsample") - 1] # type: ignore

# But can be overridden during fit
coreg_full.fit(**self.fit_params, subsample=10000, random_state=42)
assert coreg_full._meta["subsample"] == 10000
assert coreg_full.meta["subsample"] == 10000
# Check that the random state is properly set when subsampling explicitly or implicitly
assert coreg_full._meta["random_state"] == 42
assert coreg_full.meta["random_state"] == 42

# Test subsampled vertical shift correction
coreg_sub = coreg(subsample=0.1)
assert coreg_sub._meta["subsample"] == 0.1
assert coreg_sub.meta["subsample"] == 0.1

# Fit the vertical shift using 10% of the unmasked data using a fraction
coreg_sub.fit(**self.fit_params, random_state=42)
# Do the same but specify the pixel count instead.
# They are not perfectly equal (np.count_nonzero(self.mask) // 2 would be exact)
# But this would just repeat the subsample code, so that makes little sense to test.
coreg_sub = coreg(subsample=self.tba.data.size // 10)
assert coreg_sub._meta["subsample"] == self.tba.data.size // 10
assert coreg_sub.meta["subsample"] == self.tba.data.size // 10
coreg_sub.fit(**self.fit_params, random_state=42)

# Add a few performance checks
coreg_name = coreg.__name__
if coreg_name == "VerticalShift":
# Check that the estimated vertical shifts are similar
assert abs(coreg_sub._meta["vshift"] - coreg_full._meta["vshift"]) < 0.1
assert abs(coreg_sub.meta["shift_z"] - coreg_full.meta["shift_z"]) < 0.1

elif coreg_name == "NuthKaab":
# Calculate the difference in the full vs. subsampled matrices
Expand All @@ -176,7 +175,7 @@ def test_subsample(self, coreg: Callable) -> None: # type: ignore

elif coreg_name == "Tilt":
# Check that the estimated biases are similar
assert coreg_sub._meta["coefficients"] == pytest.approx(coreg_full._meta["coefficients"], rel=1e-1)
assert coreg_sub.meta["fit_params"] == pytest.approx(coreg_full.meta["fit_params"], rel=1e-1)

def test_subsample__pipeline(self) -> None:
"""Test that the subsample argument works as intended for pipelines"""
Expand All @@ -185,14 +184,14 @@ def test_subsample__pipeline(self) -> None:
pipe = coreg.VerticalShift(subsample=200) + coreg.Deramp(subsample=5000)

# Check the arguments are properly defined
assert pipe.pipeline[0]._meta["subsample"] == 200
assert pipe.pipeline[1]._meta["subsample"] == 5000
assert pipe.pipeline[0].meta["subsample"] == 200
assert pipe.pipeline[1].meta["subsample"] == 5000

# Check definition during fit
pipe = coreg.VerticalShift() + coreg.Deramp()
pipe.fit(**self.fit_params, subsample=1000)
assert pipe.pipeline[0]._meta["subsample"] == 1000
assert pipe.pipeline[1]._meta["subsample"] == 1000
assert pipe.pipeline[0].meta["subsample"] == 1000
assert pipe.pipeline[1].meta["subsample"] == 1000

def test_subsample__errors(self) -> None:
"""Check proper errors are raised when using the subsample argument"""
Expand Down Expand Up @@ -267,7 +266,7 @@ def test_coreg_raster_and_ndarray_args(self) -> None:
)

# Validate that they ended up giving the same result.
assert vshiftcorr_r._meta["vshift"] == vshiftcorr_a._meta["vshift"]
assert vshiftcorr_r.meta["shift_z"] == vshiftcorr_a.meta["shift_z"]

# De-shift dem2
dem2_r = vshiftcorr_r.apply(dem2)
Expand Down Expand Up @@ -511,19 +510,17 @@ class TestCoregPipeline:
@pytest.mark.parametrize("coreg_class", [coreg.VerticalShift, coreg.ICP, coreg.NuthKaab]) # type: ignore
def test_copy(self, coreg_class: Callable[[], Coreg]) -> None:

# Create a pipeline, add some metadata, and copy it
# Create a pipeline, add some .metadata, and copy it
pipeline = coreg_class() + coreg_class()
pipeline.pipeline[0]._meta["vshift"] = 1
pipeline.pipeline[0]._meta["shift_z"] = 1

pipeline_copy = pipeline.copy()

# Add some more metadata after copying (this should not be transferred)
pipeline._meta["resolution"] = 30
pipeline_copy.pipeline[0]._meta["offset_north_px"] = 0.5
# Add some more .metadata after copying (this should not be transferred)
pipeline_copy.pipeline[0]._meta["shift_y"] = 0.5 * 30

assert pipeline._meta != pipeline_copy._meta
assert pipeline.pipeline[0]._meta != pipeline_copy.pipeline[0]._meta
assert pipeline_copy.pipeline[0]._meta["vshift"]
assert pipeline.pipeline[0].meta != pipeline_copy.pipeline[0].meta
assert pipeline_copy.pipeline[0]._meta["shift_z"]

def test_pipeline(self) -> None:

Expand All @@ -538,8 +535,8 @@ def test_pipeline(self) -> None:
# Make a new pipeline with two vertical shift correction approaches.
pipeline2 = coreg.CoregPipeline([coreg.VerticalShift(), coreg.VerticalShift()])
# Set both "estimated" vertical shifts to be 1
pipeline2.pipeline[0]._meta["vshift"] = 1
pipeline2.pipeline[1]._meta["vshift"] = 1
pipeline2.pipeline[0].meta["shift_z"] = 1
pipeline2.pipeline[1].meta["shift_z"] = 1

# Assert that the combined vertical shift is 2
assert pipeline2.to_matrix()[2, 3] == 2.0
Expand Down Expand Up @@ -643,9 +640,9 @@ def test_pipeline_pts(self) -> None:
pipeline.fit(reference_elev=ref_points, to_be_aligned_elev=self.tba)

for part in pipeline.pipeline:
assert np.abs(part._meta["offset_east_px"]) > 0
assert np.abs(part.meta["shift_x"]) > 0

assert pipeline.pipeline[0]._meta["offset_east_px"] != pipeline.pipeline[1]._meta["offset_east_px"]
assert pipeline.pipeline[0].meta["shift_x"] != pipeline.pipeline[1].meta["shift_x"]

def test_coreg_add(self) -> None:

Expand All @@ -657,7 +654,7 @@ def test_coreg_add(self) -> None:

# Set the vertical shift attribute
for vshift_corr in (vshift1, vshift2):
vshift_corr._meta["vshift"] = vshift
vshift_corr.meta["shift_z"] = vshift

# Add the two coregs and check that the resulting vertical shift is 2* vertical shift
vshift3 = vshift1 + vshift2
Expand Down Expand Up @@ -685,21 +682,21 @@ def test_pipeline_consistency(self) -> None:
aligned_dem, _ = many_vshifts.apply(self.tba.data, transform=self.ref.transform, crs=self.ref.crs)

# The last steps should have shifts of EXACTLY zero
assert many_vshifts.pipeline[1]._meta["vshift"] == pytest.approx(0, abs=10e-5)
assert many_vshifts.pipeline[2]._meta["vshift"] == pytest.approx(0, abs=10e-5)
assert many_vshifts.pipeline[1].meta["shift_z"] == pytest.approx(0, abs=10e-5)
assert many_vshifts.pipeline[2].meta["shift_z"] == pytest.approx(0, abs=10e-5)

# Many horizontal + vertical shifts
many_nks = coreg.NuthKaab() + coreg.NuthKaab() + coreg.NuthKaab()
many_nks.fit(**self.fit_params, random_state=42)
aligned_dem, _ = many_nks.apply(self.tba.data, transform=self.ref.transform, crs=self.ref.crs)

# The last steps should have shifts of NEARLY zero
assert many_nks.pipeline[1]._meta["vshift"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[1]._meta["offset_east_px"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[1]._meta["offset_north_px"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[2]._meta["vshift"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[2]._meta["offset_east_px"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[2]._meta["offset_north_px"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[1].meta["shift_z"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[1].meta["shift_x"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[1].meta["shift_y"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[2].meta["shift_z"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[2].meta["shift_x"] == pytest.approx(0, abs=0.02)
assert many_nks.pipeline[2].meta["shift_y"] == pytest.approx(0, abs=0.02)

# Test 2: Reflectivity
# Those two pipelines should give almost the same result
Expand Down Expand Up @@ -763,7 +760,7 @@ def test_blockwise_coreg(self, pipeline: Coreg, subdivision: int) -> None:
coreg.BlockwiseCoreg(step=coreg.VerticalShift, subdivision=1) # type: ignore

# Metadata copying has been an issue. Validate that all chunks have unique ids
chunk_numbers = [m["i"] for m in blockwise._meta["step_meta"]]
chunk_numbers = [m["i"] for m in blockwise.meta["step_meta"]]
assert np.unique(chunk_numbers).shape[0] == len(chunk_numbers)

transformed_dem = blockwise.apply(self.tba)
Expand Down
Loading

0 comments on commit c30dbec

Please sign in to comment.