_agro.standardized_precipitation_evapotranspiration_index creates NaNs

[WMLKalthof]

Generic Issue

• xclim version: 0.44.0
• Python version: 3.10.12
• Operating System: Windows 10 Enterprise 22H2

Description

I am using ISIMIP chelsa-w5e5v1.0 30 arcsec resolution daily data from 1979 to 2016 (dimensions, time: 13880, y: 30, x: 57, bbox: [73.533, 18.967, 74.008, 19.217]) pr, tasmin and taxmax data to first calculate the pet (xci.potential_evapotranspiration) and the water budget (xci.agro.water_budget). No values are lost in these calculations and there are no NaNs present. However, when I then use xci.agro.standardized_precipitation_evapotranspiration_index to calculate the SPEI, two things happen:

1. Along the x axis, after a certain moment the SPEI values/pixels always become NaNs. These values do generally respond to the lowest SPEI values alongside the x-axis (in the west there is most rainfall due to orographic uplift, the eastern areas are substantially drier), however, there are lower water_budget values with other y coordinates that are not NaNs during other timesteps. Is it possible that at these pixels the algorithm is not able to fit distributions? If so, what would be a fix?
2. When I change the window below 12 periodically (corresponding to roughly 12 months) 3-5 timesteps (months) become NaNs. And with a window of 12 the first 10 timesteps become NaNs.

The settings for the function are as follows:

wb_cal = water_budget.sel(time=slice('1981-01-01', '2010-01-01'))

xci._agro.standardized_precipitation_evapotranspiration_index(wb = water_budget, wb_cal = wb_cal, freq = "MS", window = 12, dist = 'gamma', method = 'APP')

What I Received

For example, timestep 54 and 94:

And with a window of 6:

Comparison of water budget and SPEI at the same timestep:

Does anyone know what could be going wrong?

Code of Conduct

• I agree to follow this project's Code of Conduct

[huard]

Have you seen #1416 and could it be related ?

If not, I would suggest trying to use the L-Moments method to fit the gamma distribution. This is not implemented at the moment though, but should be a simple fix if you're willing to tweak xclim (I suspect it would only require adding "PWM" to the dist_and_methods dict at the top of the spi function.

[coxipi]

Hi, thanks for reporting this.

In this case #1416, the erroneous values obtained were ±np.inf, not NaN, so I'm not sure this is related. Also, the surrounding values of SPEI are not that extreme. Just to describe #1416, this is not a big computation error, essentially for very unlikely values (0.000001 or 0.999999 in the CDF), once precision limits imposed by float32 make it so that we have probabilities 0.0 or 1.0, we get ±np.inf. This precision limit is addressed in #1311 by capping the values for standardized indices to reflect this precision limit.

however, there are lower water_budget values with other y coordinates that are not NaNs during other timesteps. Is it possible that at these pixels the algorithm is not able to fit distributions?

Ok, interesting (but those lower water_budget could still translate to larger SPEI values).

The fact may still remain that you have very negative values of water_budget somewhere in the calibration for those pixels. The gamma fitting expects positive values. We have a fixed value shifting the water budget towards positive values, but this might not be enough. Could you try shifting the water_budget

wb = wb - 1.01*wb.min()

and see if the NaNs persist?

[WMLKalthof]

Have you seen #1416 and could it be related ?

If not, I would suggest trying to use the L-Moments method to fit the gamma distribution. This is not implemented at the moment though, but should be a simple fix if you're willing to tweak xclim (I suspect it would only require adding "PWM" to the dist_and_methods dict at the top of the spi function.

L-moments are invalid for this dataset it seems.

Could you try shifting the water_budget

This seems to have fixed the issue! Thank you very much for your swift assistance!

[coxipi]

Thanks to you and your detailed report, it helped troubleshooting!

Just some last remarks for future discussions about this:

1. The SPEI function in xclim already shifts values by small amount, but your case seems to reveal this was not enough. The reason for choosing a fixed amount instead of something like I just proposed e.g. : wb = wb - 1.01*wb.min() is that it can change a little the values SPEI. Say you calibrate with [1981,2010] that your whole dataset is [1981,2050], then you will have a given wb.min(). But if you extend your whole dataset to [1981,2100], then this could potentially change the wb.min() and change the shifting. So now, SPEI values in the original period [1981,2050] could be modified a little (not that much, in principle the shifting should not affect the values of SPEI, but it can still have a little influence, which is not desirable).

2. In Improve SPI performance #1311, there will be an argument in SPEI function so that the user can choose the shifting value.

3. I believe the loc parameter would allow us to deal automatically with negative values in the gamma fitting, but this may bring its share of complications. This current way of proceeding in xclim was inspired by the monocongo library, but maybe I should look more into the R package (made the researcher who proposed SPEI).