Add terrain attributes (TRI, TPI, Roughness, Rugosity), add methods of slope, aspect, hillshade and test robustness against GDAL and richDEM

[rhugonnet]

Summary of final draft

This PR adds several attributes:

• The Terrain Ruggedness Index (TRI) (also in GDALDEM),
• The Topographic Position Index (TPI) (also in GDALDEM),
• The Roughness (also in GDALDEM),
• The Fractal roughness (only in the related publication),
• The Rugosity (only in the related publication),
• The Slope, Aspect, and Hillshade can now be derived based on a refined approach (default method of GDALDEM), in addition to the method of plane fit already implemented (second method of GDALDEM, used for curvatures).

I refer to the new documentation for more detailed description and references 😉.
All the above are now tested to be exactly the same results as GDALDEM within a small tolerance (10e-3 the magnitude of the attribute). For rugosity, specific tests using simplified cases are used, and another single test is based on the example illustrated in the publication itself.

This PR also adds a wrapper around richDEM, which can be used to generate attributes available there (slope, aspect, curvatures). It is also now used to rigorously test the curvatures, that do not exist in GDALDEM, and added to the dev dependencies.

The PR also fixes some existing issues:

• There was an error in the calculation of the Planform curvature, which was just the opposite of Profile curvature, and Maximum curvature was thus also erroneous,
• The signs of Profile and Planform curvature were opposite of that of RichDEM, now corrected,
• The Hillshade differed from that of GDALDEM due to 0 being used as a no-data value, now corrected,
• The default Aspect of flat slopes in RichDEM differs from that of GDALDEM, we here use the default of GDALDEM for consistency (180° aspect).
• The default option of edge and filling methods is now "none" and "none" for consistency with other packages, and all tests have been updated consequently.

Finally, the PR improves the documentation and testing:

• All the attributes are now generated in the gallery example, and the resulting images are also displayed directly in the documentation page to avoid simply repeating ~10 times the same mini-gallery thumbnail,
• An additional multi-plot is added to the terrain attribute gallery example and used as new thumbnail,
• An additional example gallery is added on the slope and aspect methods,
• Doctests have been refined or added for all attributes,
• The reference to the publications is made for all user functions, and the exact place of the calculations (Equation number or Page number) is specified where those are computed,
• The units of all attributes is reported, except if those are unitless.
• On top of what is listed above, some tests have been added to check errors raised.

Resolves #168
Resolves #169
Resolves #170
Resolves #127
Resolves #201
Resolves #263

In terms of structure, the terrain attributes functions are based on either

• get_quadric_coefficients for 3x3 windows where calculations requires resolution,
• get_windowed_indexes for NxN windows where calculations do not require resolution,

This structure could change in time (for efficiency, potentially), but at this stage I deemed it satisfactory as performance is not a big issue.

Old comments and to-do-list (edited down)

To add on Erik's terrain features. Now we'd have everything that also exists in GDAL!
I also added references to the papers that first came up with these calculations everywhere.
For the TRI, TPI and roughness, I coded the cases where the window size can be adjusted, so not only 3x3 but any size can be used as input :).

I'm getting some differences in the tests right now (not major, but unexplained as yet). It was already the case for the slope and aspect of Erik. There exist two main algorithms for those, and we have the means of calculating the two.
Next steps:

• Include both algorithms for slope and aspect calculations,
• Make the tests against GDAL more rigorous,
• Add maximum curvature,
• Also add richDEM methods,
• Test curvature against richDEM,
• Add rugosity
• Add test for rugosity,
• Move rugosity to quadric function instead of windowed,
• Add warning for rugosity window size
• Add documentation on the new terrain attributes
• Add comment on hillshade calculation
• Add test if pixel size is not square (warning if error unclear)
• Fix tests issues due to package and geoutils updates
• Find origin of data gaps in curvatures?
• Add fractal roughness
• Test fractal roughness

[erikmannerfelt]

Wow, amazing work @rhugonnet !!

What was the small issue that fixed the GDAL difference?

Also, would all analyses be possible with your new get_windowed_indexes() or is the old function still necessary? Right now it seems like both are used, so maybe it would be good to remove the old function, if possible?

[erikmannerfelt]

Could you also add documentation entries for the new terrain attributes?

https://xdem.readthedocs.io/en/latest/terrain.html

[rhugonnet]

Yes, added in the to-do-list.
I'll do all of this when I have time again 😅

[rhugonnet]

When starting this PR again, there was a couple of failures in tests that had nothing to do with terrain.py despite the tests passing 6 months ago.
I managed to solve most of them (e.g. deprecation warnings, etc), but there is still 2 tests that are failing which are directly related to the warp_dem function and I couldn't fully grasp the source. In the meantime I added a skip decorator, and I will open a related issue so that @erikmannerfelt can check this later.

[rhugonnet]

EDIT: this is solved, see below.

Last remark of importance, especially for @erikmannerfelt
I commented the:

        # warnings.simplefilter("error")

of the tests in test_terrain.py.
I did this because they were always triggered by RuntimeWarning of NaNs, which are now present in the tests because the edge_method is "none" by default. Thus many attributes (e.g., aspect, hillshade) raise this warning when applying the functions np.arctan and others.
Not sure what we should do here to be consistent in the package? Do we silence the RuntimeWarning within the function? Or do we compute the results differently to avoid having the warning in the first place?

[rhugonnet]

Last remark of importance, especially for @erikmannerfelt I commented the:

        # warnings.simplefilter("error")

of the tests in test_terrain.py. I did this because they were always triggered by RuntimeWarning of NaNs, which are now present in the tests because the edge_method is "none" by default. Thus many attributes (e.g., aspect, hillshade) raise this warning when applying the functions np.arctan and others. Not sure what we should do here to be consistent in the package? Do we silence the RuntimeWarning within the function? Or do we compute the results differently to avoid having the warning in the first place?

This is solved, I filtered the warning as it was done for other attributes (e.g. curvatures).

[adehecq]

Phew, what a work ! Well done !! 👏
I honestly did not have the time to check each line of code, but overall it seems well documented and well tested. I'm happy to merge and add issue if they arise later.