Add a general-purpose apply_matrix() function for DEMs

[erikmannerfelt]

In #71, we discussed errors associated to DEM gridding, since ICP (and possibly more future approaches) need to convert the DEM into a point cloud. ICP currently does:

1. Convert DEM to point cloud (preparation)
2. Match the ref and tba point clouds (estimates the rigid transform)
3. Re-grid the tba DEM after transforming it (applies the rigid transform)

Step 3 currently uses scipy.interpolate.griddata which is extremely slow! A faster and more reliable approach would be fantastic to have. Something with the approximate signature:

def apply_matrix(dem: np.ndarray, transform: rio.transform.Affine, matrix: np.ndarray) -> np.ndarray:
    # do matrix magic
    return transformed_dem

One idea I had was to describe the 3D rigid transform as a 2D similarity transform and a linear vertical correction:

A 3D rigid transform from above (from 2D) would include horizontal translation, horizontal rotation, horizontal scaling (due to rotation-induced perspective differences), and elevation (bias + rotation) differences. As far as I understand, the first 3 parts could be described as a 2D similarity transform, and the final as a deramp.

The advantage of such an approach is that the DEM never has to be converted to a point cloud and regridded; it just needs to be moved and twisted a bit. I guess the most difficult part of this, however, would be to mathematically do the conversion...

This general-purpose function is integral for #79 to implement.

I am also all-ears to other alternatives.

[erikmannerfelt]

In my sleep, I came up with a hacky solution that should work for now:

1. Transform a set of arbitrary points using the 3D rigid transform.
2. Estimate a 2D (horizontal) SimilarityTransform using skimage.transform.estimate_transform(orig_points2D, trans_points2D, "similarity")
3. Estimate a (vertical) deramp using the trusty xdem tools
4. Apply the SimilarityTransform using skimage and then apply the deramp.

The great advantage of this is that the input DEM is never converted into a point cloud, and thus doesn't need regridding, but is just warped (skimage) and value corrected (deramp).

Of course, the best would be to estimate the SimilarityTransform and Deramp mathematically instead of empirically deriving them, but that's beyond my math knowledge!

[erikmannerfelt]

Added in #87