Bug with raster-point coregistration

Hi, I was about to co register some IceSat2 data with a reference raster using the NuthKaab Method. When I used the apply method of the co-registration object it seems to shift the points in the opposite direction. Syntax I used looks like this:

    NuthKaab=xdem.coreg.NuthKaab()
    coreg_result=NuthKaab.fit(ref_dem, points, z_name ='z', verbose=True)
    points2=NuthKaab.apply(points)

Unfortunately, test_coreg_example_shift in test_affine .py does not capture this error, so edited it a little to make the bug reproducible:

    def test_coreg_point_shift(self, shift_px, z_shift=2, coreg_class=coreg.NuthKaab, verbose=False, subsample=5000):
        """
        For comparison of coreg algorithms:
        Shift a ref_dem on purpose, e.g. shift_px = (1,1), and then applying coreg to shift it back.
        """
        res = self.ref.res[0]

        # shift DEM by shift_px
        shifted_ref = self.ref.copy()
        shifted_ref.shift(shift_px[0] * res, shift_px[1] * res, inplace=True)

        shifted_points = shifted_ref.to_pointcloud(
            subsample=subsample, force_pixel_offset="center", random_state=42
        ).ds
        points_on_ref = self.ref.to_pointcloud(
            subsample=subsample, force_pixel_offset="center", random_state=42
        ).ds

        shifted_points["E"] = shifted_points.geometry.x
        shifted_points["N"] = shifted_points.geometry.y
        shifted_points.rename(columns={"b1": "z"}, inplace=True)
        points_on_ref.rename(columns={"b1": "z"}, inplace=True)
        shifted_points["z"]=shifted_points["z"]+z_shift

        kwargs = {} if coreg_class.__name__ != "GradientDescending" else {"subsample": subsample}

        coreg_obj = coreg_class(**kwargs)

        coreg_obj.fit(self.ref, shifted_points, verbose=verbose, random_state=42)
        corrected_points = coreg_obj.apply(shifted_points)

        if ((corrected_points.geometry.x-points_on_ref.geometry.x)[0] == pytest.approx(0, abs=0.1)) and ((
                corrected_points.geometry.y-points_on_ref.geometry.y)[0] == pytest.approx(0, abs=0.1)) and ((
                    corrected_points['z']-points_on_ref['z'])[0]== pytest.approx(0, abs=0.1)):
            return
        east_diff = (corrected_points.geometry.x-points_on_ref.geometry.x)[0]
        north_diff = (corrected_points.geometry.y-points_on_ref.geometry.y)[0]
        print(f"Diffs are too big. east: {east_diff:.2f} m, north: {north_diff:.2f} m")

I traced the cause of this behavior to coreg/affine.py in the method : NuthKaab._fit_rst_pts. In line 1014 to 1015 there is the following code:

self._meta["shift_x"] = offset_east * resolution if ref == "point" else -offset_east 
self._meta["shift_y"] = offset_north * resolution if ref == "point" else -offset_north

I don't see a reason why the shift for points is stored in pixels. From what I see, the if condition here is unnecessary and deleting it will let the test I created run smoothly. Also it is a little confusing with the signs. They actually should be there so it makes sense. But, since the signs are opposite in the _apply_pts function they don't need to be taken care off in _fit_rst_pts. It was a bit hard to get my head around this and it would be probably better to have the apply method consistent for points and rasters as well as for height and lateral shift.

This could also concern other co registration methods, where the exact same lines are found.

I hope I didn't get anything wrong with the syntax and the behavior was intended as I assumed.

Thanks for opening this issue @JonasLiebsch!

I did have some doubts about the points methods for the affine classes, given how poor the tests were when this was added a year ago :sweat_smile:, which I why I opened this: https://github.com/GlacioHack/xdem/issues/485 and tried to not document it in the documentation before solving it.

I have been tackling it here in the past weeks: https://github.com/GlacioHack/xdem/pull/530, where I've deleted all _apply_ functions for Affine classes, which means the code will always use consistently the same apply_matrix function that is tested quite thoroughly :). I still need to finalize a couple other aspects, and add the same robust test suite as for point-raster methods in BiasCorr, and then it should all be good.

Good catch for the offset error! This might have been a mistake from me in the recent re-work of the Coreg metadata (https://github.com/GlacioHack/xdem/pull/526), which would have come up if the point tests for Affine classes were better :sweat_smile:

I'm leaving for long holidays at the end of the week, so I might not have the time and only pick that up again in August. If so, I'll try to take an hour or two to at least fix this issue specifically :wink:

GlacioHack / xdem

Bug with raster-point coregistration #546