Improve mapping interpolation with k>1 when NaN values are present

[ortk95]

Currently, we use RectBivariateSpline, so have to replace NaN pixels with some value before interpolating (currently just set to 0). Even though we later mask any projected NaN pixels (with _should_propagate_nan_to_map), this can produce edge effects near NaN pixels when using e.g. cubic splines.

It would make more sense therefore to either:

• Replace NaNs with a more appropriate constant value before interpolating - e.g. the median value of the image rather than 0. This would be the easiest to implement, and have negligible performance impact.
• Replace NaNs with locally derived value before interpolating - e.g. do nanmean of 3x3 footprint around each pixel. This would cost more performance, but provide better results for variable data.
• Use a different interpolation method that doesn't require a regular grid - probably using bisplrep or griddata.

[ortk95]

Could do locally derived value fixing with something along the lines of

def clean(img: np.ndarray) -> np.ndarray:
    nans = np.isnan(img)
    cleaned = img.copy()
    cleaned[nans] = np.nanmedian(img)
    to_fix = nans & ~scipy.ndimage.uniform_filter(nans, size=3)
    for i, j in np.argwhere(to_fix):
        cleaned[i, j] = np.nanmean(img[max(i - 1, 0) : i + 2, max(j - 1, 0) : j + 2])
    return cleaned

[ortk95]

griddata clearly performs best at keeping original variation with no influence from NaNs, so think we should use that as the default method. Can keep RectBivariateSpline as an option for e.g. if smoothing is wanted.

[ortk95]

From futher testing, griddata doesn't seem as great:

• It can become prohibitively slow for larger (~1000px) datasets, as it needs to construct a triangulation.
• The results can be unstable and dependent on the system (see failing tests only on GitHub's macOS workers). It seems like this is probably again due to the triangulation, which doesn't necessarily have a unique solution, and can therefore be affected by floating point errors.

I think this probably therefore rules out griddata, as e.g. we don't want unstable mapping results to be interpreted as changes in the target etc.

Scipy's RegularGridInterpolator seems to perform really nicely (and quickly) when paired with our custom _replace_nans_with_interpolated_values, so this is probably the way to go.