Reported discrepancy between LS Factor in SDR vs SAGA

[phargogh]

I've heard this now from several people including Dionee on the natcap forums and now from Rafa as well: there is a significant difference between the calculated LS factor in InVEST SDR when compared with SAGA's LS factor. Rafa sent along his inputs to the SAGA module in QGIS (download them here, Stanford auth required) for further investigation.

[phargogh]

Rafa noted today that:

There is no absolute deadline for it. We are waiting on it for a paper but it's not super urgent.

[phargogh]

Doug and Jesse had a coincidentally very relevant exchange on slack this morning, where Jesse thought that slope was in radians when it turns out that pygeoprocessing produces slope in degrees. The SDR UG also assumes that slope is in radians. It could be that this LS factor discrepancy is simply due to a unit conversion issue.

[phargogh]

Nope, that's not it. _calculate_ls_factor does the correct conversion (arctan(percent_slope/100)) for percent slope to slope in radians.

[phargogh]

Ok, so I've been looking into this for a little while now and here's a summary of what I've found so far:

It's worth noting that I'm using SAGA version 7.3.0 in QGIS 3.22.4-Białowieża, while the source code I'm referencing is a development version on master, specifically: saga-gis/saga-gis@c2e0847.

SAGA's LS Factor doesn't handle plateaus

SAGA's LS factor calculation relies on a flow accumulation method, Get_Flow() that only considers a pixel downstream if and only if its elevation is strictly less than the current pixel's elevation. Thus, this method has no support for plateaus which is problematic because our pitfilling algorithm produces lots of plateaus.

The solution for this is to use the SAGA Wang & Liu pitfilling algorithm to ensure suitable hydrological connectivity. When using Wang&Liu on the DEM, the InVEST LS Factor and the SAGA LS Factor are actually on the same order. Not identical (like ~30-50% off), but at least they're roughly the same order of magnitude.

SAGA's LS Factor Uses Strict Desmet & Govers 1996; InVEST Does Not

The critical difference is in the $m$ parameter.

$$ m_{SAGA} = \frac{\beta}{1 + \beta}$$

Given:

• $\beta = \frac{\sin\theta / 0.0896}{3\sin\theta^{0.8} + 0.56}$
• $\theta$ is the slope of the current pixel for which we are calculating the LS factor in radians

Contrast that with the InVEST $m$ parameter.

$$ \begin{align*} m_{InVEST} &= \left\{\begin{array}{lr} 0.2, & \text{where } \theta \leq 1\% \\\ 0.3, & \text{where } 1\% < \theta \leq 3.5\% \\\ 0.4, & \text{where } 3.5\% < \theta \leq 5\% \\\ 0.5, & \text{where } 5\% < \theta \leq 9\% \\\ \frac{\beta}{1 + \beta}, & \text{where } \theta > 9\% \end{array}\right\} \\\ \\\ \beta &= \frac{\sin\theta / 0.0896}{3\sin\theta^{0.8} + 0.56} \end{align*} $$

This accounts for some of the difference between the two LS methods.

Saga Uses On-Pixel Aspect; InVEST Does Not

The SAGA source code uses the on-pixel aspect calculated from their Get_Gradient() function (example call), which is the method described in Zevenbergen & Thorne 1987 and repeated in Desmet & Govers 1996 and is then adjusted into an aspect-length factor $x$ given the aspect in radians, $\theta$, according to Desmet & Govers 1996:

$$ x = |sin(\theta)| + |cos(\theta)| $$

InVEST by contrast infers aspect from the weighted average of flow contributions to each neighbor according to the Multiple Flow Direction algorithm.

While pixel values have reasonable values in the range $ [1, \sqrt{2}] $ in both cases, the SAGA approach has more spatially contiguous outputs:

Contributing Area Differs between SAGA and InVEST

The SAGA LS Factor operation offers multiple methods for calculating the upstream contributing area, including:

1. Contour length simply as cell size
2. Contour length dependent on aspect

InVEST does neither of these, instead using the literal upstream contributing area and not treating it as the accumulated flow length as SAGA does.

So, is there a problem?

If you are using a Wang&Liu-conditioned DEM, these other differences are minor.

However, the aspect issue seems like something we may want to correct. I'm not sure about the others; I'll check with Rafa about this.

[phargogh]

I just met with Rafa and we're going to continue to do some exploratory development and data analysis that will probably end up with some modifications to the InVEST LS Factor function. For now, there are a few things that need to be done:

• Document the need for using a Wang & Liu-conditioned DEM when comparing LS factor outputs from SAGA with LS Factor outputs from InVEST. Otherwise it just won't be an apples-to-apples comparison because the flow accumulation will be very off.
• Rafa is very in support of using true aspect in the LS factor instead of the weighted MFD aspect.
• Rafa is also in support of using the trivial upstream length in InVEST instead of the upstream area as we do now.
• I will write a script to condition the DEM with SAGA's Wang & Liu and use the resulting DEM to run the SAGA LS factor on each of the upstream area methods as well as InVEST using each of these upstream area methods.

[phargogh]

I forgot to update this. We're trying out a few different methods for updating the LS factor on the Costa Rica inputs that Kelly used for a recent project and Rafa will do the data analysis to figure out what's the best solution for InVEST's SDR from there.

[phargogh]

As a reminder to myself, I'm working on this in https://github.com/phargogh/invest-ls-factor-vs-saga

[phargogh]

We're waiting on some real-world data to run this on. The original hope was to use the Costa Rica layers, but maybe we should use something else?

[phargogh]

I just chatted with @schmittrjp about this on slack and here's his response:

I think let's put it in the latest release. That work on CR has stalled as Becky has moved on, but from what I have seen (in CR, and Peru) the updated LS factor is the way to go.

So, @dcdenu4 this is a bit of a bigger change. Once the change is ready, should the release wait until 3.14?

@dcdenu4 Would it be OK with you if I added an ADR for this? (related to #1079)

[dcdenu4]

@phargogh, 3.14 probably makes sense. Given how much this update effects results we might want to have a forum post as well with some visuals outlining the change.

I think an ADR would be awesome. I know I had done some research and prep for an updated proposal related to #1079 . I'll see if I can dig any of that out... Actually I think I was writing up an ADR for the Habitat Quality change. I'll see where that ended up.

[phargogh]

I have just created a PR for this, #1289 and I am just waiting on final confirmation from Rafa about the changes we're hoping to make before submitting the PRs for this and for an ADR discussing the change.

[phargogh]

The ADR discussing this change is at #1288 and can be previewed here: https://github.com/phargogh/invest/blob/task/1079-915-adr-for-ls-factor-change/doc/decision-records/ADR-0001-Update-SDR-LS-Factor.md

[phargogh]

Word from @schmittrjp is:

I'd say this change has now been well tested in a number of settings and outperformed the current LS factor. So let's include it in the next release!