Python segmentation fault on Windows/Linux

[groutr]

I am working on large DEM sized (11999, 20000) with a dtype of float32 and experiencing a segfault on windows and linux.

rdDEM = richdem.rdarray(DEM, no_data=raster.nodata)
rdDEM.geotransform = transform.to_gdal()
richdem.FillDepressions(rdDEM, epsilon=args.epsilon_fill, in_place=True)
richdem.ResolveFlats(rdDEM, in_place=True)

print("Flow Direction computation...")
flow = richdem.FlowProportions(rdDEM, method='D8')    # <--- segfault happens here
print("done")

This was the backtrace I was able to obtain from the coredump.

(gdb) bt
#0  0x00002b68d5e9ffa0 in void richdem::FM_OCallaghan<(richdem::Topology)1, float>(richdem::Array2D<float> const&, richdem::Array3D<float>&) ()
   from xxxxxxx/lib/python3.11/site-packages/_richdem.cpython-311-x86_64-linux-gnu.so
#1  0x00002b68d5f890c7 in void pybind11::cpp_function::initialize<void (*&)(richdem::Array2D<float> const&, richdem::Array3D<float>&), void, richdem::Array2D<float> const&, richdem::Array3D<float>&, pybind11::name, pybind11::scope, pybind11::sibling, char [5]>(void (*&)(richdem::Array2D<float> const&, richdem::Array3D<float>&), void (*)(richdem::Array2D<float> const&, richdem::Array3D<float>&), pybind11::name const&, pybind11::scope const&, pybind11::sibling const&, char const (&) [5])::{lambda(pybind11::detail::function_call&)#3}::_FUN(pybind11::detail::function_call&) () from /scratch2/NCEPDEV/ohd/Ryan.Grout/mambaforge/envs/rivershapes/lib/python3.11/site-packages/_richdem.cpython-311-x86_64-linux-gnu.so
#2  0x00002b68d5ef3952 in pybind11::cpp_function::dispatcher(_object*, _object*, _object*) ()
      from xxxxxxx/lib/python3.11/site-packages/_richdem.cpython-311-x86_64-linux-gnu.so
#3  0x000055db28859c56 in cfunction_call ()
#4  0x000055db28839bd3 in _PyObject_MakeTpCall.localalias ()
#5  0x000055db28846463 in _PyEval_EvalFrameDefault ()
#6  0x000055db288fccbd in _PyEval_Vector ()
#7  0x000055db288fc34f in PyEval_EvalCode ()
#8  0x000055db2891b12a in run_eval_code_obj ()
#9  0x000055db28917433 in run_mod ()
#10 0x000055db2892c070 in pyrun_file ()
#11 0x000055db2892ba1e in _PyRun_SimpleFileObject.localalias ()
#12 0x000055db2892b634 in _PyRun_AnyFileObject.localalias ()
#13 0x000055db28925baf in Py_RunMain.localalias ()
#14 0x000055db288eb007 in Py_BytesMain ()
#15 0x00002b68cb588555 in __libc_start_main () from /lib64/libc.so.6
#16 0x000055db288eaead in _start ()

Interestingly, richdem doesn't segfault when using a slightly smaller DEM. I tried 11500x20000 and it ran successfully. I think the system has plenty of memory available because watching memory usage under windows, at the time of the crash, there is still a healthy amount of free memory.

EDIT:
This is running Python 3.11 and richdem 2.3.0.

[groutr]

I was able to build my own version of richdem and produce a coredump with a generated test case.
The backtrace is pointing to this particular line.
https://github.com/r-barnes/richdem/blob/c5a183ada30fbd1b8587121e0d313b4cda9390c0/include/richdem/flowmet/OCallaghan1984.hpp#L26C3-L26C29

#0  richdem::Array3D<float>::setAll (val=-1, this=0x558b2871c5c0) at xxxxxxx/richdem-dev/x86_64-conda-linux-gnu/include/c++/12.3.0/bits/unique_ptr.h:191
#1  richdem::FM_OCallaghan<(richdem::Topology)0, float> (elevations=..., props=...) at xxxxxxx/richdem-dev/include/richdem/flowmet/OCallaghan1984.hpp:26

[groutr]

Test case:

import numpy as np
import richdem

k = np.random.rand(11999, 20000).astype('float32')
rd = richdem.rdarray(k, no_data=-999)
rd.geotransform = (1.0, 0, 0, 0, 1.0, 0)

richdem.FlowProportions(rd, method='D8')

[groutr]

After playing with this some more, I'm wondering if maybe the culprit is in the python bindings. I'm able to process my DEM with the standalone executables that richdem builds.