Use libtopotoolbox to fill sinks

[wkearn]

I just merged TopoToolbox/libtopotoolbox#14, which implements sink filling for DEMs.

Let's figure out what we need to do to expose this to users through the Python library.

The function fillsinks has the signature

void fillsinks(float *output, float *dem, ptrdiff_t nrows, ptrdiff_t ncols);

where output and dem are both arrays of floats with size nrows * ncols. I think implementing the bindings to call this function with Numpy arrays or a future GridObject class (#6) should be fairly straightforward.

It seems more important at this point to make sure we can build and access the library from Python on different platforms because then we can start more confidently adding functionality to the library. If we need to make some changes to the build process for libtopotoolbox, for example, it will be easier to make those now before it gets too complicated.

Let's come up with a plan here and open some more issues either here or in TopoToolbox/libtopotoolbox to keep track of more specific things we need to do to get this working.

@Teschl: what do you think?

[wschwanghart]

Very cool. I didn't test anything so far, but looking at the code of fillsinks.c, I was wondering, whether the program can handle nans.

[wkearn]

That is a good question. It does not do anything specific to handle NaNs, which means the filled DEM probably won't be correct. I'd have to think about what would actually happen though. What should be the correct behavior when a DEM has NaNs in it?

[Teschl]

I have been experimenting around a bit, and just to confirm: all functions will be built into the one libtopotoolbox.so file?

One potential problem I encountered was that ptrdiff_t is not supported in ctypes. I've read that it's possible to simply use integers, but I haven't been able to properly test this yet because I'm still trying to figure out how to correctly build the package while including the library.

[wkearn]

Yes, the idea is that the whole library will be provided as a single shared library file.

ptrdiff_t is likely the same size as ssize_t, which does seem to be available in ctypes. I am not particularly tied to using ptrdiff_t for the array sizes, if it would make it significantly easier to use something different.

[Teschl]

~~I think I would prefer using something other than ptrdiff_t. It doesn't crash when using ssize_t, size_t, or just long/int, but it just keeps running for minutes. It does complete if I set nrows/ncols to 10, for example, so it's hard to know exactly what is going wrong.~~

[Teschl]

Correction: The reason for the unexpected behavior was the presence of NaNs in the TIF I was using. Without NaNs, it works just like the libraries I tested on before. Also, ssize_t seems to work just fine.

[wkearn]

Yes, I'm realizing that the NaNs are a problem at the moment. It runs for a long time because the sink filling never converges when there are NaNs, and I handled the maximum iteration cutoff incorrectly (cf TopoToolbox/libtopotoolbox#20). I think I have some fixes for it, though.

[wkearn]

With #14, we can now run

import topotoolbox as tt3
dem = tt3.GridObject("dem.tif")
demfs = dem.fillsinks()

from Python to fill sinks using the morphological reconstruction algorithm implemented in libtopotoolbox.