Preserve nodata in masked_array or Raster subclass during Coreg operations

[rhugonnet]

EDIT with reproducible example (of NaNs created by Coreg but not masked), while both reference_raster and to_be_aligned_raster are fine:

import xdem
import geoutils as gu
import numpy as np

xdem.examples.download_longyearbyen_examples()

reference_raster = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
to_be_aligned_raster = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_tba_dem"))
glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
inlier_mask = ~glacier_mask.create_mask(reference_raster)

nuth_kaab = xdem.coreg.NuthKaab()
nuth_kaab.fit(reference_raster.data, to_be_aligned_raster.data,
              inlier_mask=inlier_mask, transform=reference_raster.transform)
aligned_raster = nuth_kaab.apply(to_be_aligned_raster)
diff = reference_raster - aligned_raster

print(np.min(diff.data))
print(np.ma.min(diff.data))
print(np.ma.min(aligned_raster.data))
print(np.ma.min(reference_raster.data))
print(np.ma.min(to_be_aligned_raster.data))

nan
nan
nan
8.052481
8.376112

For opening issue: we don't want to repeat this operation after every Coreg operation to preserve nodata along calculations:

aligned_raster.data[~np.isfinite(aligned_raster)] = reference_raster.nodata
diff = reference_raster - \
           gu.Raster.from_array(aligned_raster.data, transform=reference_raster.transform, crs=reference_raster.crs, nodata=reference_raster.nodata)

Originally posted by @rhugonnet in #224 (comment)

[rhugonnet]

@adehecq I think this is closed by GlacioHack/geoutils#257, right?

[rhugonnet]

Actually, it's not as it happens within the Coreg classes. Though similar issue as #178 and #145

[adehecq]

So I'm not sure I understand. Can it be simplified now or not?
I believe you don't need the call to gu.Raster.from_array anymore if you provide a Raster directly to a Coreg method following #175.
Does the arithmetic overloading not work?

[rhugonnet]

I have updated the top comment with a reproducible example.

@adehecq For keeping a written trace on this: arithmetic overloading works! The issue is about the propagation of NaNs during Coreg operations, and in the mask of the masked_array.
I think the following lines in Coreg.apply() might be the issue as they simply copy the old mask instead of making a new one, so it's an issue for both masked_arrays and RasterType:

# If the DEM was a masked_array, copy the mask to the new DEM
        if hasattr(dem, "mask"):
            applied_dem = np.ma.masked_array(applied_dem, mask=dem.mask)  # type: ignore
        # If the DEM was a Raster with a mask, copy the mask to the new DEM
        elif hasattr(dem, "data") and hasattr(dem.data, "mask"):
            applied_dem = np.ma.masked_array(applied_dem, mask=dem.data.mask)  # type: ignore

        # If the input was a Raster, return a Raster as well.
        if isinstance(dem, gu.Raster):
            return dem.from_array(applied_dem, transform, dem.crs, nodata=dem.nodata)

As we just discussed, this is also linked to GlacioHack/geoutils#231 and other various nodata issues
in geoutils

[adehecq]

I see. I get a similar issue with nodata handling (#239). I'll try to look into that and solve it this afternoon if I can.
But it seems like apply_matrix is not returning NaNs where it should.