Proposal: terrain attributes based on 5x5 polynomial fit

[trchudley]

Posted here following discussion with @rhugonnet

Currently, xdem calculates terrain coefficients over a 3x3 window following Zevenbergen and Thorne (1987) and Horn (1981). This is suitable approach for the vast majority of scenarios, but there are potential alternatives which have advantages in specific situations.

An alternative approach by Florinsky (2009) calculates derivatives by fitting a third-order partial polynomial to a 5x5 window. This larger window can be better suited to high-resolution photogrammetry data as it leads to a local denoising effect.

I have implemented this method in a numba-friendly way in my own work for ArcticDEM/REMA strips and would be happy to help porting this functionality to xdem if desired - if only because in the medium term I would like to move my coregistration functions to use xdem in the backend!

P.S. Elsewhere in Florinsky's work, there are also a number of further curvature attributes that might be desirable to port into xdem and are easy to implement. This includes horizontal, vertical, and gaussian curvature as well as unsphericity.

[rhugonnet]

Amazing, your implementations are much shorter, probably also more computationally efficient than ours (I wonder if our "if" loops slow down the process or not at all, I don't have a good enough understanding of numba)! 😅

For this type of kernel/array size relation, I remember reading a blogpost showing that convolution was typically faster than numba, which is why we have this ongoing: #486.
But we need to do some more testing before actually implementing it. It could be the perfect occassion to compare with your numba implementation too, then move forward from here 🙂