NuthKaab shifted grids in fit_func contain interpolated nodata

[rhugonnet]

Sorry for opening so many issues on the Nuth and Kaab in the same day, @erikmannerfelt! It means we are testing (sometimes struggling) thoroughly 😉
The Nuth and Kaab shift function based on scipy.interpolate.RectBivariateSpline has a wrong NaN propagation that keeps interpolated NaNs.
Reproducible example:

import numpy as np
dem = np.arange(12.).reshape((4, 3))
dem[1, 1]= np.nan
print(dem.astype(int))

[[ 0. 1. 2.]
[ 3. nan 5.]
[ 6. 7. 8.]
[ 9. 10. 11.]]

import scipy.interpolate

# Make index grids for the east and north dimensions
east_grid = np.arange(dem.shape[1])
north_grid = np.arange(dem.shape[0])

# Make a function to estimate the aligned DEM (used to construct an offset DEM)
elevation_function = scipy.interpolate.RectBivariateSpline(
    x=north_grid, y=east_grid, z=np.where(np.isnan(dem), -9999, dem), kx=1, ky=1
)

# Make a function to estimate nodata gaps in the aligned DEM (used to fix the estimated offset DEM)
# Use spline degree 1, as higher degrees will create instabilities around 1 and mess up the nodata mask
nodata_function = scipy.interpolate.RectBivariateSpline(
    x=north_grid, y=east_grid, z=np.isnan(dem), kx=1, ky=1
)

offset_east = 0.5
offset_north = 0.5

origin_elevation = elevation_function(y=east_grid, x=north_grid)
print(origin_elevation.astype(int))

[[ 0 1 2]
[ 3 -9999 5]
[ 6 7 8]
[ 9 10 11]]

new_elevation = elevation_function(y=east_grid + offset_east, x=north_grid - offset_north)

# Set NaNs where NaNs were in the original data
new_nans = nodata_function(y=east_grid + offset_east, x=north_grid - offset_north)
new_elevation[new_nans >= 1] = np.nan

print(new_nans)

[[0. 0. 0. ]
[0.25 0.25 0. ]
[0.25 0.25 0. ]
[0. 0. 0. ]]

print(new_elevation.astype(int))

[[ 0 1 2]
[-2498 -2497 3]
[-2495 -2494 6]
[ 8 9 9]]

The fix should be quite simple, replacing:

new_elevation[new_nans >= 1] = np.nan

by

new_elevation[new_nans > 0] = np.nan

[rhugonnet]

Based on similar arguments (interpolating something with -9999 or -32XXX always yields something aberrant), I would apply the same threshold 0 everywhere, including:

• coreg.apply_matrix(), currently:

tr_nan_mask = skimage.transform.warp(
      nan_mask.astype("uint8"),
      inds,
      order=resampling_order,
      mode="constant",
      cval=1,
      preserve_range=True
  ) > 0.1  # Due to different interpolation approaches, everything above 0.1 is assumed to be 1 (True)

• coreg.warp_dem() (if not deprecated), currently:

new_mask = skimage.transform.warp(
    image=dem_mask,
    inverse_map=np.moveaxis(new_indices, 2, 0),
    output_shape=dem_arr.shape,
    cval=False
) > 0.25

I might have missed other occurences!

[adehecq]

Damn, this coregistration really seem to cause a lot of troubles with NaNs!
This is somewhat related to PR #197 and the associated issue. In there @erikmannerfelt changed the threshold for gaps to 0.1 but I can't remember why exactly...