Bindings and xtensor-compatible API for core features of LSDTT

[benbovy]

I'm interested in the possibility of ultimately reusing core features of LSDTT directly from other C++ applications or from Python, i.e, bypassing the drivers, I/O files, etc. It's great to have drivers and Python tools built on top of it, but it has also some overhead and limited flexibility.

Besides all the cmd-line tools provided by LSDTT, it would be great to have an API that would look as simple as the Python snippet below:

import numpy as np
import lsdtt

dem = np.load(...)   # load a DEM as a numpy array
floodplains = lsdtt.extract_floodplains(dem)    # directly calls the C++ code

Maybe I'm missing something and this already exists? If not, I think the xtensor library would help to achieve this.

xtensor is a relatively new C++ array library that is designed very much in the style of numpy. I use it for fastscapelib and I'm pretty happy with it!

It is still a young project but it has the following assets:

• The developer team is very active, responsive and cares about the design.
• There is a growing ecosystem around it (xtensor-blas, xtensor-fftw, xsimd...).
• Once we have an xtensor-compatible C++ API, it is easy to also provide bindings in Python/R/Julia (thanks to xtensor-python, xtensor-r and xtensor-julia).
• xtensor's array data model is extensible and can wrap pretty much any array container. I've written converters for xtensor arrays <-> richdem's Array2d objects and it was easy. I guess it wouldn't be to too hard writing converters for the TNT array class(es) or the LSDArray2D class.
• LSDTT and fastscapelib would interoperate seamlessly :-)

Any thoughts? Note that this could be maintained in a 3rd-party repository (e.g., something like xtensor-lsdtt), but this would require careful handling of the dependencies.

I'll try writing a simple wrapper function for some lsdtt feature as a proof of concept when I have some time.

[fclubb]

I think this would be a really great idea, and would potentially help uptake of the code a lot as well as providing more flexibility for users. Let me know if you want to chat about it at some point @benbovy, I think some of the drivers that do simple slope and curvature maps plus channel extraction would be a good place to start.

[simon-m-mudd]

I fully support this. It doesn't exist already and I've not tried xtensor but sounds promising. Boris has done some experimenting with numba. Is there an advantage of xtensor over numba?

[dvalters]

@bgailleton and myself have chatted a lot about how to go about writing Python bindings to topotools. We reasoned it would be a lot of work to get a working Python API, but I haven't come across xtensor before (we mainly discussed cython and/or numba). I think it is a very cool idea though.

[bgailleton]

Numba was not flexible enough to write bindings, but with cython it worked pretty well. Basically we (@dvalters helped me a lot) wrote a binding in c++ that translate the tnt 2d arrays into vectors (cython understand vectors) and compiled it as a shared library that we can then call from python Exactly the way you are describing. I have never used xtensors but it can worth a try? At the moment I just ported the polyfir function to cython and still need to find how to compile a DLL for Windows so if xtensors is more straightforward it can be good.

[benbovy]

I don't really know if one day it will be possible/easy to call C++ compiled code from numba jitted functions, but right now if you want to use numba then I think that you have to re-implement the features in Python, and possibly maintain two implementations (one in C++ and one in Python). Just-in-time compilation with numba has also pros and cons (more flexible but may take some time when importing the library to compile the functions that you want to use).

Using something like xtensor, you only need to write some wrappers around the existing LSDTT C++ code.

Actually, xtensor is not the only alternative for writing Python bindings. As @bgailleton says, Cython could do a good job here too. You could also use pybind11 directly, which has some support for numpy arrays.

Unlike Cython or pybind11, the xtensor ecosystem already provides tools to create bindings in other languages as well. You may find it useful if you want to later call LSDTT routines from R or Julia (I haven't tried this yet, though). xtensor may also help if you want to provide a top-level API in C++. I find it pretty cool to be able to write code in C++ that looks like numpy. Depending on how xtensor will be adopted in the future, it will be (hopefully) easy to inter-operate with other C++ libraries (but who knows).