Define consistent nodata propagation for interp_points and allow to return interpolator

[rhugonnet]

Most of the line changes in this PR are due to re-structuring, moving out functions of the Raster class (see #539, in preparation for #383) into their own module, here only those linked to interp_points.

This PR also adds a nodata propagation scheme (as scipy.interpolate.interpn does not support nodata propagation natively except for methods "nearest" and "linear") and allows to return an interpolator that contains this nodata propagation (to optimize speed of iterative methods calling the same interpolator several times, such as DEM coregistration).

It turned out that interpolation of map_coordinates with splines of order 3/5 does not match "cubic" or "quintic" in interpn, only order 0/1 do match for "nearest" and "linear", so all other are dropped. The nodata propagation and accuracy of results for values near the edge of NaNs is now tested thoroughly! 😄

Resolves #559
Starts addressing #539 for functions related to interp_points

TO-DO:

• Add tests on return_interpolator,
• Add tests that gets exact interpolation value, then creates artificial NaNs, and compares to the interpolation near NaN edges to estimate the accuracy,
• Adjust nodata propagation tests for "nearest" and "linear" (by enforcing the same propagation between map_coordinates and interpn)?