Use euclidean distance to propagate NaNs in interp_points?

[rhugonnet]

Currently we use a series of binary_dilation to propagate NaNs to a certain pixel distance in interp_points, but this seems to match the L1 norm (abs(dX) + abs(dY)) and not a classic euclidean distance.
Should we change this scheme to an euclidean distance instead? (same calculation as Raster.proximity())
Would this better match the reason for the nodata propagation? (=size of interpolation kernels for linear/cubic/quintic methods)

See original mention below:

Originally posted by @rhugonnet in #560 (comment)