How to calculate the SPEI index only with daily precipitation and temperature datasets?

[CGL5230]

Setup Information

• Xclim version:0.46.0
• Python version:3.11.0
• Operating System:Linux

Context

Hi,all! I want to ask how to calculate the SPEI index only with daily precipitation and temperature datasets?

I find a function is xclim.indices.standardized_precipitation_evapotranspiration_index in the xclim. But there is no water budget dataset. I find another funtion is xclim.indices.water_budget in the xclim which need more parameters.

I read the notebook of climate indices. It can calculate the SPEI, only with pr and tas dataset. But, it seems like only fit to the statio scale. That means the Need for Loop Statements and the cost for the time and source.

Is there any way can meet my need in xclim? Thanks a lot for your help :)

Code of Conduct

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

[coxipi]

Hi!

notebooks and loops and performance

I'm not sure which notebook you are referring to, I can't find any notebook using water_budget or SPEI functions, what's the name of the notebook?

Is it also related to what you say about "loop" and "station scale"? You don't need to do loops over location points, this will be handled by SPEI function. Just like if you have an input dataset with (time,lat,lon), you can compute da.max(dim="time"), and this will return you an array (lat,lon), you don't need to loop over specific location points. For SPEI, inputting your water budget will return you an array with original dimensions (time,lat,lon) (although the size of time dimension may be reduced, say if you input a daily dataset but require a monthly SPEI)

That being said, the SPEI function can get quite resource heavy, so as you scale things up, you may need to be careful with your resource management. You might want to look at dask which allows to separate your dataset in chunks and control how much RAM is allowed to be used: https://docs.xarray.dev/en/stable/user-guide/dask.html

Let me know if I misunderstood your questions.

using water_budget

To use only tas and pr to get water budget, you need to use the method "MB05". The various methods are documented in potential_evapotranspiration function. This method computes evapotranspiration using only tas and the extraterrestrial solar radiation obtained from the latitudes in your dataset, so latitudes need to be in tas or supplemented to your water_budget call. Use this code:

wb= xclim.indices.water_budget(pr=pr, tas = tas, method="MB05")

Then you can use the wb as input for the SPEI function xclim.indices.standardized_precipitation_evapotranspiration_index(wb,...)

comments about water_budget errors

I will say that water_budget is somewhat complicated to use with its many optional arguments. You CAN use only tas, but in this specific case, you NEED to specify the correct method "MB05". I think that the function implied potential_evapotranspiration could give a clearer error message, say if you give some input variables but they're not compatible with the method selected (or the default one).

Currently, with the failing code

wb = xclim.indices.water_budget(pr=pr, tas = tas)

we get the error

line 1378, in potential_evapotranspiration
    tasmin = convert_units_to(tasmin, "degF")
line 427, in convert_units_to
    raise NotImplementedError(f"Source of type `{type(source)}` is not supported.")

which is apparently the error you get for inputting None in convert_units_to. By itself, I don't think this message is super effective in letting me know "oops, you inputed None where it's unexpected". More importantly, a higher level error message at the level potential_evapotranspiration could help users learn about their error in using evap methods.

[CGL5230]

Thanks @coxipi ！

I try to calculate the wb followed your suggestion as wb_AWI_CM_ssp126= xclim.indices.water_budget(pr=pr_AWI_CM_ssp126.pr, tas =tas_AWI_CM_ssp126.tas, method="MB05") .
Then, I use SPEI_AWI_CM_ssp126=xclim.indices.standardized_precipitation_evapotranspiration_index(wb=wb_AWI_CM_ssp126,freq='MS',window=6) to calculate SPEI-6.

I got some warnings:
'RuntimeWarning: Mean of empty slice.' ' m = x_pos.mean()' and 'RuntimeWarning: invalid value encountered in scalar divide
ret = ret.dtype.type(ret / rcount)'

I got the SPEI-6 finally,like below. I think I followed the standard steps to calculate SPEI index. The result will be correct.