TopoToolbox 3 Expands Quantitative Geomorphology Analysis Across Python and MATLAB Platforms

Researchers have released TopoToolbox 3, restructuring the established terrain analysis platform around an independent computational core to bring quantitative geomorphology workflows to Python alongside MATLAB. Published in Earth Surface Dynamics by William S. Kearney and an international research team, the update expands a software suite that has supported landscape analysis across terrestrial geomorphology, planetary geology, glaciology, and hydrology for fifteen years. Earlier versions operated exclusively within proprietary MATLAB environments, requiring specialized toolboxes that restricted access and complicated integration with outside software. The new release addresses these limitations by introducing a multi-layer framework designed to broaden accessibility, facilitate integration with external simulation tools, and adopt formal research software engineering standards under the FAIR principles. Essential computational algorithms have been reorganized into a shared library, ensuring reproducible geomorphic calculations across programming environments while preserving full backward compatibility for existing users.

At the foundation of this redesigned architecture is libtopotoolbox, an independent library written in C that implements core computational routines through a data-oriented design. Rather than managing complex object state or internal data ownership, libtopotoolbox operates strictly on contiguous memory arrays passed directly by host languages. Language-specific binding layers—pybind11 for Python and MEX for MATLAB—marshal array data between user sessions and the underlying C routines. The high-level user interface preserves the object-oriented structure established in previous editions, organized around Grid objects representing digital elevation models, Flow objects describing flow directions, and Stream objects representing channel networks. In MATLAB, these retain familiar class names like GRIDobj, FLOWobj, and STREAMobj, while the new Python package exposes corresponding GridObject, FlowObject, and StreamObject classes. Both implementations represent flow networks as topologically sorted edge lists, ensuring that all incoming edges to any network node occur before outgoing edges, which permits rapid linear scans to evaluate upstream contributing areas and drainage metrics.

Deploying scientific algorithms across distinct languages created significant quality assurance demands, particularly surrounding the oracle problem, where true analytical solutions for complex terrain operations on real topography are rarely known in advance. To verify code without relying on subjective inspection, the developers constructed an automated testing framework comprising snapshot characterization testing, property-based testing, and metamorphic evaluation. Snapshot tests capture outputs produced by previous MATLAB versions across sets of real and synthetic digital elevation models, storing them as reference baselines. Automated test suites in libtopotoolbox and Python download these reference pairs, recompute results using their own implementations, and confirm that outputs match within narrow floating-point limits. In addition, property-based testing verifies mathematical invariants over randomly generated elevation surfaces without needing known global answers. For instance, tests automatically confirm that sink-filling functions eliminate all closed depressions and that drainage basin delineation assigns identical basin identifiers to any pixel and its downstream neighbor.

The testing pipeline also relies heavily on metamorphic testing, which identifies subtle defects by examining whether expected relationships hold between outputs when input data undergo systematic transformations. Because specific operations should yield predictable outcomes—such as maintaining identical slope magnitudes when an elevation grid is rotated by ninety degrees—violations uncover hidden implementation flaws. In one practical demonstration, metamorphic tests detected an inconsistency in libtopotoolbox’s steepest-descent flow routing across different memory layouts. MATLAB formats multidimensional grids in column-major order, whereas NumPy in Python defaults to row-major structures. Although the core library processed arrays correctly in memory, the directional scan examining neighboring cells searched clockwise for column-major matrices and counterclockwise for row-major matrices. When encountering equal elevations among adjacent cells, this discrepancy caused the algorithm to resolve ties in opposing directions, altering flow paths and producing differing flow accumulation rasters. Adding an ordering parameter synchronized neighborhood scanning directions across layouts, ensuring identical flow paths regardless of memory organization.

To illustrate typical scientific workflows across both language interfaces, the researchers implemented a comparative analysis computing Chi maps for the Rhine and Danube drainage basins in southwestern Germany and northern Switzerland. Rooted in steady-state stream power models of river incision, the Chi metric integrates the inverse of upstream contributing drainage area along channel reaches, transforming longitudinal river profiles into linear plots that expose transient imbalances driven by tectonic uplift or river capture. The workflow loaded Copernicus Global Digital Elevation Model data through OpenTopography queries, projected the raster into an ETRS89-LAEA coordinate system at ninety-meter resolution, and removed sub-zero elevations to define base level. Flow routing followed D8 steepest descent with least-cost topographic carving to traverse sinks, after which stream networks were initiated at a ten-square-kilometer contributing threshold. Isolating the two primary river networks through connected-component filtering enabled calculation of along-channel Chi values, which were visualized spatially using platform-specific plotting utilities.

Evaluating the resulting Chi profiles revealed minor differences between the two environments, with the Python implementation underestimating Chi along the main trunks of the Rhine and Danube relative to the MATLAB output. Diagnostic testing proved that this discrepancy did not stem from calculation routines in libtopotoolbox, which executes identical flow routing, accumulation, and integration algorithms across both implementations. Instead, the deviation originated in external reprojection libraries used during spatial preprocessing. Reprojecting raw elevation rasters from geographic coordinates to the projected grid relies on MATLAB’s Mapping Toolbox in one environment and Python’s rasterio library in the other. Although both pipelines nominally applied bilinear resampling to identical spatial resolutions, differing internal interpolation assumptions produced slightly different interpolated elevation matrices. When the team eliminated this external variable by feeding pre-projected elevation datasets directly into both packages, the resulting Chi profiles matched perfectly, highlighting how external preprocessing dependencies can influence computational reproducibility across scientific platforms.

Beyond static river network metrics, TopoToolbox 3 incorporates hydrodynamic modeling capabilities through the integration of GraphFlood, a numerical method developed to estimate basin-scale flow physics efficiently. While one-dimensional channel representations are effective for regional analyses, they do not capture two-dimensional variables such as inundation width, water depth, flow velocity, and bed shear stress, whereas solving full two-dimensional shallow-water equations requires substantial computational resources. GraphFlood balances physical realism and computational tractability by combining graph theory with simplified shallow-water equations to resolve two-dimensional steady-state flow over complex landscapes. To showcase its performance within TopoToolbox 3, the authors simulated flow in the Saison watershed of the Western French Pyrenees using a twenty-five-meter resolution digital elevation model from the French National Institute of Geographic and Forest Information. Effective precipitation rates of 10, 50, and 100 millimeters per hour were prescribed across the catchment, allowing researchers to evaluate water depth and extract shear stress, hydraulic slope, and channel width across a five-hundred-meter swath around the main river.

GraphFlood’s incorporation highlights the modular utility of libtopotoolbox, which executes the hydrodynamic core while language interfaces handle user interaction and output management. In Python, GraphFlood is encapsulated within a GFObject class that configures run parameters, orchestrates input arrays, and outputs results directly as GridObject rasters. In MATLAB, the model runs through functional calls that exchange data with existing spatial matrices. Because the hydrodynamic model takes standard Grid objects as inputs and returns simulation outputs as Grid objects, users can immediately pass results into downstream TopoToolbox functions for spatial filtering, swath analysis, or visualization. Embedding such specialized physical process models directly within the shared C library allows research teams to distribute complex numerical tools broadly across scientific communities without forcing users to compile external software or write custom file converters.

The Python binding layer also enables bidirectional coupling with external modeling frameworks, which the authors demonstrated by interfacing TopoToolbox 3 with Landlab, an open-source Python environment for simulating Earth surface dynamics. The team constructed a landscape evolution experiment combining Landlab’s stream power river incision and linear hillslope diffusion components with a custom flow routing module, TT3FlowRouter. This component substitutes Landlab’s native flow accumulation tool with TopoToolbox’s D8 routing and morphological depression-carving algorithms. Within the simulation, TT3FlowRouter converts Landlab raster fields into GridObject instances, computes upstream drainage areas, and translates TopoToolbox’s flow network into Landlab’s directional link format at every time step. After running a simulated domain to steady state over two million years, the uplift rate was doubled from one to two millimeters per year, triggering transient knickpoints that migrated upstream through the river network. TopoToolbox was then applied post-simulation to identify knickpoint locations and plot their upstream progression on hillshaded digital elevation models and Chi profiles.

Computational benchmarking during the Landlab integration revealed performance advantages as well as inter-model communication bottlenecks. Executing the landscape evolution model with TopoToolbox’s flow router required 24.1 seconds on an eight-core workstation, compared to 44.7 seconds when using Landlab’s FlowAccumulator paired with its depression-routing module. While Landlab ran faster when depression finding was disabled entirely (11.5 seconds), TopoToolbox currently resolves depressions by default during flow routing. Detailed profiling indicated that translating TopoToolbox’s topologically sorted edge list into Landlab’s network links consumed twice as much computation time as flow routing itself, underscoring the overhead of transforming network topologies across software boundaries. Furthermore, differences in depression-carving algorithms between the packages led to noticeable local elevation discrepancies. On unconditioned random initial topography, divergent flow routing caused elevation differences up to 200 meters on individual cells, though carving an initial minimum-slope channel reduced the final mean absolute elevation difference across the entire domain to 0.8 meters.

By adopting open-source distribution under the GNU General Public License v3.0, persistent Zenodo DOIs, and automated GitHub continuous integration, TopoToolbox 3 aligns computational geomorphology with the FAIR principles for research software. Active development centers on reaching full feature parity between Python and MATLAB packages, including porting swath profiling and drainage divide analysis tools, while progressing dedicated R packages and QGIS desktop plugins. To address memory constraints when evaluating massive high-resolution digital elevation models, developers are also exploring hardware acceleration frameworks such as PyTorch, Taichi Lang, and TensorFlow to enable multi-processor and GPU computing. The authors caution that automated verification tests confirm algorithmic compliance with software specifications but do not replace scientific validation against empirical field observations. By lowering barriers between programming languages and establishing transparent software development practices, TopoToolbox 3 provides an adaptable platform that invites community-driven contributions to advance quantitative research into Earth surface processes.

Subject of Research: Open-source software architecture and cross-platform terrain analysis in quantitative geomorphology

Article Title: TopoToolbox 3, a laboratory for quantitative geomorphology

Article References: Kearney, W. S., Arnau, G., Bringezu, T., Dietze, M., Gailleton, B., Lamprecht, A.-L., Lenz, K., Ott, R., Scherler, D., Steer, P., & Schwanghart, W. (2026). TopoToolbox 3, a laboratory for quantitative geomorphology. Earth Surface Dynamics, 14(6), 865-885. https://doi.org/10.5194/esurf-14-865-2026

Image Credits: AI Generated

DOI: 10.5194/esurf-14-865-2026

Keywords: quantitative geomorphology, digital elevation models, TopoToolbox, libtopotoolbox, flow routing, stream power equation, chi analysis, metamorphic testing, property-based testing, GraphFlood, Landlab, landscape evolution modeling, FAIR principles, drainage basin delineation, knickpoint migration, pytopotoolbox, hydrodynamic modeling