Refactor monthly quickflow function #1109

emlys · 2022-11-03T21:03:15Z

Description

Fixes #1318
This PR refactors the monthly quickflow function with the goal of making it clear how the different cases and edge cases are handled. I reorganized the function around 3 main cases, one of which has a couple edge cases. I tried to make the code line up as closely as possible with the quickflow equation in the user's guide. I hope that the big comment block clarifies what's going on.

The functional changes from #1318 are:

When the s_i / a_im ratio is greater than 100 (an arbitrary value I chose as per discussion with Lisa and Rafa), set QF to 0
When QF is negative, set it to 0

And a few trivial changes:

remove trailing zeros left over from python 2 era
use the same TARGET_NODATA variable for all output nodatas
remove 2 unused parameters of the monthly quickflow function

Checklist

Updated HISTORY.rst (if these changes are user-facing)
Updated the user's guide (if needed)
Tested the affected models' UIs (if relevant)

emlys · 2022-11-15T23:22:58Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

    def qfi_sum_op(*qf_values):
        """Sum the monthly qfis."""
-        qf_sum = numpy.zeros(qf_values[0].shape)
-        valid_mask = ~utils.array_equals_nodata(qf_values[0], qf_nodata)


We were assuming that all 12 monthly quickflow rasters would have nodata in the same areas, but that's not necessarily true.

Why would that not be true? If a pixel does have no quickflow in a month, I'd expect it to be 0. If its outside of the model domain it would be always nan?

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

emlys · 2022-11-15T23:35:53Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

-        exp_result[nonzero_e1_mask] = numpy.exp(
-            (0.8 * valid_si[nonzero_e1_mask]) / a_im[nonzero_e1_mask] +
-            numpy.log(E1[nonzero_e1_mask]))


The old implementation took advantage of exponent math to avoid overflow: By rewriting exp(0.8 * s_i / a_im) * E1(s_i / a_im) as the equivalent exp(0.8 * s_i / a_im + log(E1(s_i / a_im))), the result of exp wouldn't get so large. But then, there's another edge case where E1 = 0 and log(E1) = infinity. That was handled with the nonzero_e1_mask above.

I replaced that with handling the overflow warning, below. It's closer to the original equation, and IMO it's clearer what's going on.

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

dcdenu4

Hey @emlys, thanks for looking into this and being thorough. I like your approach in breaking out the components and the documentation that accompanies it. I think the handling of the OverflowError is works well, but I do wonder if there's a way to constrain the problem further up somehow. I'll be interested to hear what Rafa says.

Oh, I didn't write this as a comment, but is there a reason you chose to "ignore" the error vs surround it in a try / except block?

Thanks!

dcdenu4 · 2022-11-17T15:13:40Z

HISTORY.rst

@@ -45,6 +45,9 @@ Unreleased Changes
      now reprojected to the ``lulc_cur_path`` raster. This fixes a bug where
      rasters with a different SRS would appear to not intersect the
      ``lulc_cur_path`` even if they did. (https://github.com/natcap/invest/issues/1093)
+* Seasonal Water Yield
+    * Fixed a bug where monthy quickflow nodata pixels were not being passed
+      on to the total quickflow raster (`#1105 <https://github.com/natcap/invest/issues/1105>`_)


Since this is user facing, maybe an additional note of: This could result in negative values on the edges

dcdenu4 · 2022-11-17T17:46:29Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

+        #     Solution: Catch the overflow warning, and set the term to 0
+        #     anywhere that overflow happens.


Are there any numerical constraints that make sense for s_i, a_im, or s_i / a_im that would allow us to mask this out proactively? I'd be curious from the science side of what these terms represent when these boundaries are hit.

We're setting to 0 because if the exp has overflowed, then we know the E1 term is 0?

Values of s_i / a_im that make exp overflow aren't common, but are totally possible with reasonable input data, for example:

CN = 30 P = 8 (millimeters of rain per month) n = 10 (number of rain events per month) s_i = 1000/CN - 10 = 23.33 a_im = P/(n * 25.4) = 8/(10 * 25.4) = 0.031 s_i / a_im = 23.33 / 0.031 = 752 exp(752) = 3.8 * 10^326, which overflows float64.

as far as numerical constraints, exp(x) overflows when x is about 709.78. I'm not sure if this might vary on different systems.

We're setting to 0 because if the exp has overflowed, then we know the E1 term is 0?

yes, the E1 will be extremely close to 0 and "underflow" to 0. for example:

>>> scipy.special.exp1(200) 6.885226106307636e-90 >>> scipy.special.exp1(400) 4.776013586420972e-177 >>> scipy.special.exp1(700) 1.406518766234033e-307 >>> scipy.special.exp1(800) 0.0

@emlys thanks for making this example. It makes numerically sense, but let's look at it from a process perspective. basically we have no precip (0.8 mm per rain is basically no rain) meeting a very permeable soil (very low curve number). As I see from the above, if we just use the threshold value of exp(709) then QF would be 0, right? Probably we can check that somehow bottom up? Calculate for each pixel what values of x would result in an overflow and set quickflow to 0 there? There might e a more elegant way. Let's think about it after thanksgiving.

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

dcdenu4 · 2022-11-17T17:58:39Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

+
+        # case 1: there is no precipitation where both p_im and n_m are
+        # defined and equal to zero.
+        case_1_mask = ~precip_mask


I might be mistaken or misinterpreting but case_1_mask here could include nodata pixels, which would not be "defined".

precip_mask ~precip_mask X V 0 0 1 0 1 0 1 X = NoData V 0 0 1 0 0 0 1 1 V = Value > 0 V V X 1 1 0 0 0 1 0 = Value == 0

because precip_mask is defined above as:

# precip mask: both p_im and n_m are defined and greater than 0 precip_mask = valid_p_mask & valid_n_mask & (p_im > 0) & (n_m > 0)

it could include pixels that have nodata in the s_i or stream arrays, but not in the precip or n_events arrays.

That's because you don't need to know s_i or stream to calculate QF = 0 when P = 0. I think you could argue that it should be nodata anyway, but this is consistent with the past implementation.

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

dcdenu4 · 2022-11-17T18:15:42Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

+        term_result[subcase_a_mask] = 0
+
+        # edge case 3b: set the whole term to 0 when exp() overflows
+        subcase_b_mask = numpy.isinf(exp_result)


Interesting, so ignoring the OverflowError will leave exp_result with an infinity value?

Yes:

>>> numpy.exp(800) <stdin>:1: RuntimeWarning: overflow encountered in exp inf

To your question above

is there a reason you chose to "ignore" the error vs surround it in a try / except block?

Because it's a warning not an error, and a RuntimeWarning could be caused by other things whereas the numpy context manager lets you specifically ignore overflows

Oh interesting! I was testing math.exp(800) and getting:

>>> math.exp(900) Traceback (most recent call last): File "<stdin>", line 1, in <module> OverflowError: math range error

That makes sense.

dcdenu4 · 2022-11-17T18:32:08Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

+        case_3_mask = valid_mask & precip_mask & ~stream_mask
+
+        term_result = numpy.full(
+            qf_im[case_3_mask].shape, TARGET_NODATA, dtype=numpy.float32)


Is qf_im[case_3_mask].shape == qf_im.shape? Is there a reason to mask here?

Masking gets you the shape of where the mask is true, i.e. a 1-d array whose length is sum(mask).

>>> a = numpy.array([True, False, True, False]) >>> b = numpy.array([1, 1, 1, 1]) >>> b.shape (4,) >>> b[a].shape (2,)

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

Co-authored-by: Doug <[email protected]>

…tcap#1318

emlys · 2023-06-23T00:08:27Z

@dcdenu4 I've updated this with the latest from our conversation with @lmandle and @schmittrjp last week!

dcdenu4

Thanks @emlys,

I had a few comments / suggestions.

I was also wondering if we talked about updating the Users Guide at all to better reflect what the model was doing under some of these instances? I can't remember if during our discussion with Rafa or Lisa that came up.

dcdenu4 · 2023-06-27T13:15:12Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

+        valid_mask = (
+          valid_p_mask &
+          valid_n_mask &
+          ~utils.array_equals_nodata(stream, stream_nodata) &


It looks like this valid_mask rework is missing the stream_mask. This might work itself out below, but just noting it as I work through.

Yeah, the stream mask is separate

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

dcdenu4 · 2023-06-27T13:56:48Z

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

+        )
+
+        # case 3d: set any negative values to 0
+        qf_im[qf_im < 0] = 0


Is the nodata value we're using here positive or is it -1? I'm seeing TARGET_NODATA set as -1, so maybe we should have a nodata mask here?

Thanks for catching that!

dcdenu4 · 2023-06-27T14:11:58Z

tests/test_seasonal_water_yield_regression.py

+            expected_quickflow_array, atol=1e-5)
+
+    def test_monthly_quickflow_large_si_aim_ratio(self):
+        """Test `_calculate_monthly_quick_flow` with undefined nodata values"""


Update docstring here?

dcdenu4

Thanks @emlys, let's roll this in!

emlys added 7 commits November 3, 2022 14:01

simplify some things in SWY quick flow function

16c47eb

remove unneeded divide-by-zero mask

bc53960

better handle edge cases in SWY quick flow function

bd0dbde

use target nodata consistently, remove trailing .0s

8590203

reorganize quickflow sum function

487150d

organize quickflow code by three cases

1ed6df6

finish organizing quick flow cases and write explanatory comment

4ce6789

emlys commented Nov 15, 2022

View reviewed changes

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py Show resolved Hide resolved

emlys commented Nov 15, 2022

View reviewed changes

src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py Outdated Show resolved Hide resolved

emlys requested review from dcdenu4 and schmittrjp November 16, 2022 00:02

emlys self-assigned this Nov 16, 2022

add history note

3608df2

emlys changed the title ~~SWY quick flow function~~ Fix quickflow masking bug and refactor quickflow function Nov 16, 2022

Merge branch 'main' into bugfix/1105

dd7292f

emlys marked this pull request as ready for review November 16, 2022 00:09

dcdenu4 requested changes Nov 17, 2022

View reviewed changes

Update src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py

89199e3

Co-authored-by: Doug <[email protected]>

dcdenu4 added the on hold There's a reason we're not working on this yet label May 8, 2023

emlys changed the title ~~Fix quickflow masking bug and refactor quickflow function~~ Refactor quickflow function Jun 2, 2023

remove changes directly related to natcap#1105

16154b3

emlys changed the base branch from main to feature/output-spec June 2, 2023 17:59

emlys changed the base branch from feature/output-spec to main June 2, 2023 17:59

Merge branch 'main' into bugfix/1105

59aa1fd

emlys removed the on hold There's a reason we're not working on this yet label Jun 21, 2023

emlys added 3 commits June 22, 2023 16:51

clean up refactor, handle case where QF is negative, and add tests na…

6f219d9

…tcap#1318

Merge branch 'main' into bugfix/1105

9f154c3

set very small expected values to 0 in tests

b7388c1

add history note natcap#1318

248d1f1

emlys requested a review from dcdenu4 June 23, 2023 00:08

emlys changed the title ~~Refactor quickflow function~~ Refactor monthly quickflow function Jun 23, 2023

Merge branch 'main' into bugfix/1105

07325e2

dcdenu4 requested changes Jun 27, 2023

View reviewed changes

emlys and others added 2 commits June 27, 2023 13:17

Merge branch 'main' into bugfix/1105

ac0d73b

fix test docstring; qf_im nodata masking natcap#1105

399f7b0

emlys requested a review from dcdenu4 July 19, 2023 18:54

dcdenu4 approved these changes Jul 31, 2023

View reviewed changes

dcdenu4 merged commit 8045806 into natcap:main Jul 31, 2023

dcdenu4 mentioned this pull request Jul 31, 2023

SWY quickflow update natcap/invest.users-guide#128

Closed

emlys deleted the bugfix/1105 branch October 3, 2024 23:00

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Refactor monthly quickflow function #1109

Refactor monthly quickflow function #1109

emlys commented Nov 3, 2022 •

edited

Loading

emlys Nov 15, 2022

schmittrjp Nov 18, 2022

emlys Nov 15, 2022

dcdenu4 left a comment

dcdenu4 Nov 17, 2022

dcdenu4 Nov 17, 2022

emlys Nov 17, 2022

emlys Nov 17, 2022

schmittrjp Nov 18, 2022

dcdenu4 Nov 17, 2022

emlys Nov 17, 2022

dcdenu4 Nov 17, 2022

emlys Nov 17, 2022

dcdenu4 Nov 17, 2022

dcdenu4 Nov 17, 2022

emlys Nov 17, 2022

emlys commented Jun 23, 2023

dcdenu4 left a comment

dcdenu4 Jun 27, 2023

emlys Jul 19, 2023

dcdenu4 Jun 27, 2023

emlys Jul 19, 2023

dcdenu4 Jun 27, 2023

dcdenu4 left a comment

		# Solution: Catch the overflow warning, and set the term to 0
		# anywhere that overflow happens.

Refactor monthly quickflow function #1109

Refactor monthly quickflow function #1109

Conversation

emlys commented Nov 3, 2022 • edited Loading

Description

Checklist

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dcdenu4 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

emlys commented Jun 23, 2023

dcdenu4 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dcdenu4 left a comment

Choose a reason for hiding this comment

emlys commented Nov 3, 2022 •

edited

Loading