Repository navigation
Discrepancies in flow routing between libtopotoolbox and TopoToolbox v2 #122
Description
I noted in #72
There are differences between the two on flat regions that have to do with how the auxiliary topography is treated. TT2 itself computes different flow directions over flats with its MEX implementation and the default MATLAB implementation. The flow directions of the libtopotoolbox implementation seem more similar to the MATLAB implementation than to the MEX version. There is an unused function in src/flow_routing.c, compute_flowdirection_TT2, which attempts to more directly copy the MEX implementation. However, using compute_flowdirection_TT2 still matches the MATLAB implementation better than it matches the MEX implementation. I am open to exploring these discrepancies more, but I would like to have something that we can start building from.
The discrepancies I noticed when testing #72 seemed quite minor and limited to the flats, but these differences get amplified when computing the drainage basins from the flow routing information using #121.
Here are two maps of drainage basins on the Big Tujunga DEM included with TopoToolbox v2. The left one is computed with
DEM = GRIDobj("DEMdata/srtm_bigtujunga30m_utm11.tif");
FD = FLOWobj(DEM);
L = drainagebasins(FD);while the right one is computed using libtopotoolbox's flow routing.
It may be hard to tell, but among other things, there are twice as many basins in the libtopotoolbox version than in the TopoToolbox v2 version.
Most or maybe all of those extra basins are on the boundary, which makes me think that in addition to having problems routing over flats, we might also not be identifying outlets on the boundary of the DEM correctly. I didn't notice those earlier because I was primarily looking at maps of flow accumulation, and these small basins have very little flow accumulation.
I would like to figure out what is going on here before merging #121.
Activity
wkearn commentedon Sep 3, 2024
Maybe this is a better figure. Top is TopoToolbox v2 and bottom is libtopotoolbox.
Implement snapshot tests
Add a snapshot test for the auxiliary topography
Change CarveMinVal to a double
Change CarveMinVal to a double (#155)
wkearn commentedon Jan 6, 2025
Some experimentation has convinced me that this is not due to the implementation of
drainagebasins, so I merged #121.The flow routing problems may at least partially be related to numerical errors that accumulate in the gray-weighted distance transform results. This is the difference between the auxiliary topography computed using the MATLAB implementation (
graydist) and libtopotoolbox (gwdt) on the Tibet DEM:These differences are small (the maximum error is around 0.09), and I still need to confirm whether it affects the flow routing in flat regions. The MATLAB implementation converts the costs to double precision floating point numbers before running
graydist, and this seems to be the source of the differences in the resulting auxiliary topography. When you rungraydistwith single precision floats, there are still some numerical differences between the libtopotoolbox and the MATLAB implementations, but they remain small and do not really accumulate over the flat region. I am not yet sure if this means we need to use double precision distances in ourgwdtimplementation or not.Add Drainage basin delineation section to the Getting Started guide
Add Drainage basin delineation section to the Getting Started guide (#…
wkearn commentedon Jan 8, 2025
The numerical problems seem not to be the largest source of these flow routing discrepancies. Instead there are two related problems, which only affect flow routing over flat areas:
Flow routing over constant slope auxiliary topography
How should the flow routing algorithm distinguish between pixels with identical values for the auxiliary topography? Say your auxiliary topography looks like this
which pixel should the central pixel (value 3) be routed towards?
steepestneighbor_mexwill choose the northeastern neighbor, because it only looks at the minimum value of the auxiliary topography not the gradient -- it is not scaled by the inverse square root of two. If more than one neighboring pixel has the minimum value of the auxiliary topography, it will choose the first pixel clockwise from north.FLOWobjwill choose the southeastern neighbor, because it first sorts the pixels by the auxiliary topography then by the real topography and then chooses the /largest linear index/ of a neighboring pixel, which is always to the south and east of the central pixel in column major ordering.These three different options result in stream networks that look like
respectively.
Sills on the border
The other problem is how sills that lie on the edge of the DEM are treated. This is the reason you see the weird patterns in the drainage basins in this image from the Tibet test case:
Imagine a DEM that looks like this:
All of the pixels on the eastern boundary are sills:
A sill on a boundary can always be an outlet in any of the three implementations, however, the different implementations end up doing different things:
steepestneighbor_mexboth do the following: a sill is not a flat, so we look at the DEM values rather than the auxiliary topography. Each of our sills is surrounded by neighbors that are either higher or at the same elevation as it, so it cannot drain anywhere and must be an outlet. Each sill becomes its own outlet. Because of which neighbor the two different algorithms choose, libtopotoolbox creates five drainage basins:and
steepestneighbor_mexcreates three:with the two south east pixels disconnected from the flow network. These drain independently of each other, so one might also label them as single-pixel drainage basins, but the current flow routing with
steepestneighbor_mexdoes not pick them up.In the Tibet DEM, the flat area in the east ends up with drainage basins that look like this for the three implementations:
Potential solutions
The flow directions in libtopotoolbox can easily be made to match
steepestneighbor_mexeasily: we need to start looking at neighbors from the north and not scale by the inverse square root of two for the auxiliary topography. Both could also be made to match the native MATLAB implementation by scanning through neighbors counterclockwise from the southeast. I'm not sure there is a good scientific argument to prefer one or another of these flow directions. I like the libtopotoolbox implementation because it treats the auxiliary topography just like the real topography. You switch to using the auxiliary topography when you hit a flat.The boundary conditions are a bigger problem. Our current approach to boundaries in flow routing is very ad hoc, and this gets us into the situation above, where different algorithms provide very different answers because of how they handle the boundary. It would be nice to write down exactly how we intend to treat different boundaries. This particular problem visible in the Tibet DEM occurs when sills are on the boundary, but I am not confident that we handle all other boundary conditions consistently either.
wkearn commentedon Jan 8, 2025
This situation does not actually require the boundary pixels to be labeled as sills. Instead it occurs when neighboring boundary pixels share an elevation but have no lower elevation pixel nearby. This happens in the Big Tujunga DEM. There is a region on the southern boundary that looks like
which both libtopotoolbox and
steepestneighbor_mexassign to drainage basins asEach pixel at elevation 489 is made its own outlet. The native MATLAB implementation labels everything as one drainage basin because all the elevation 489 pixels flow into each other to the east.
So the boundary pixels that are a problem are those that have no lower neighbor. Libtopotoolbox and
steepestneighbor_mexboth label every boundary pixel that has no lower neighbor as its own outlet. The native MATLAB implementation instead routes flow over the boundary pixels with no lower neighbors by choosing the neighbor with the highest linear index.Implement libtopotoolbox-based flow routing