SPEI : can't figure it out #1589
-
Setup Information
ContextHello, I have been trying to calculate the SPEI for a long time unsuccessfully, because every time I tried to fix an error I got a new one . I've been using agERA5 dataset, rearranged from GRIB to netCDF, adding the "pev" variable, that is the potential evapotranspiration calculated with the Penman-Monteith equation. The netCDF file contains several variables with daily values (for february 2023, as an example only), including those necessary for computing the SPEI: in particular the "ppn" (that is the daily total precipitation [mm]) and the "pev" (namely the daily potential evapotranspiration [mm]) in order to compute the daily hydroclimatic balance [mm] (precipitation minus evapotraspiration) named "bic_daily". Below the first (simplest) version of the python code i wrote: I got this error: KeyError: 'units' (same error if I change "indices" with "atmos"). The error des not depends on the parameters inside the SPEI function. I decided to try with some changes. Below the new python code: If I run: I decided to try again with some other changes. Below the last python code: And I got this: AttributeError: 'DataArray' object has no attribute 'time' I do not understand if the problem comes from the data, because I have not found examples of code for the SPEI computation Giulia Code of Conduct
|
Beta Was this translation helpful? Give feedback.
Replies: 10 comments
-
|
Sorry to hear you rushed using SPEI. As you perform an operation on an arithmetic operation on your datasets in xarray, ( # a dataset with units
pr.attrs["units"]
>>> 'kg m-2 s-1'
# arithmetic op kills units
out = 0.5*pr
out.attrs["units"]
>>> KeyError: 'units'
# specify that ops should not change units
with xr.set_options(keep_attrs=True):
out = (0.5*pr)
out.attrs["units"]
>>> 'kg m-2 s-1'If you wrap your computation of As for the I guess if you solve the problem in the first step, you won't be doing this change |
Beta Was this translation helpful? Give feedback.
-
|
Hi @coxipi, thank for your answer. I did what you said, but I get "KeyError: 'units'". I think I have to set some attributes manually because in the xarray.datarray I use for the SPEI computation there are no attribute at all. I don't understand which kind of attributes I need to set. |
Beta Was this translation helpful? Give feedback.
-
|
Hi! Can you show the code you use to get this error? It seems you might have invoked |
Beta Was this translation helpful? Give feedback.
-
|
Yes it's not defined in my initial dataset. How can I define it? And how can I define every attributes I need to compute the SPEI? My code is exactly the code I wrote above, plus the code you advise me to add. |
Beta Was this translation helpful? Give feedback.
-
|
water budget is Now, suppose the units are agERA5_202302['bic_daily'] = agERA5_202302['ppn'] - agERA5_202302['pev']
bic_daily_foramonth = agERA5_202302['bic_daily']
bic_daily_foramonth.attrs["units"] = "kg m-2 s-1"You could also define units for agERA5_202302['ppn'].attrs["units"] = "kg m-2 s-1"
agERA5_202302['pev'].attrs["units"] = "kg m-2 s-1"If you find out that the implicit units aren't the same, you can use our conversion tools to change units of one datasets agERA5_202302['ppn'].attrs["units"] = "kg m-2 s-1"
agERA5_202302['pev'].attrs["units"] = "kg mm-2 d-1"
agERA5_202302['pev'] = xclim.core.units.convert_units_to(agERA5_202302['pev'], agERA5_202302['ppn']) # 'pev' is now in "kg m-2 s-1"Cheers |
Beta Was this translation helpful? Give feedback.
-
|
Hi @coxipi and thank you very much for your time. Here is my new code: agERA5_202302 = xr.open_dataset('202302_agERA5.nc') And here is my new error: Traceback (most recent call last): During handling of the above exception, another exception occurred: Traceback (most recent call last): |
Beta Was this translation helpful? Give feedback.
-
|
It's impossible I can't compute the SPEI. Is anyone here who can tell me why I'm getting a problem with the reindex_time function?. It seems to be expecting a variable or dimension named 'day' in the dataset, but it's not finding it, even if I checked the dataset and everithing seems correct to me. Another code: import xarray as xr
import xclim.indices as indices
from xclim.indicators import atmos
agERA5_202302 = xr.open_dataset('202302_agERA5.nc')
agERA5_202302['bic_daily'] = agERA5_202302['ppn'] - agERA5_202302['pev']
bic_daily_foramonth = agERA5_202302['bic_daily']
bic_daily_foramonth.attrs["units"] = "kg m-2 s-1"
bic_daily_foramonth.attrs['frequency'] = 'D' # Assuming it's a daily frequency
spei_feb = indices.standardized_precipitation_evapotranspiration_index(wb = bic_daily_foramonth, freq = None, window = 2, dist =
'gamma', method = 'APP')Alaways getting the same error: My time coordinates: Thank you in advance and sorry to keep writing |
Beta Was this translation helpful? Give feedback.
-
|
Hi @GiuliaIsNotAvailable, indeed you are hitting a bug within xclim, we're sorry about that! The "day" that raises the error should be "dayofyear". We'll try to fix it as soon as possible. Sadly, that means that without editing xclim's code itself, it is currently impossible to compute a daily SPI/SPEI... |
Beta Was this translation helpful? Give feedback.
-
|
Hi @aulemahal. I shall wait with confidence a code update. Thank you very much. |
Beta Was this translation helpful? Give feedback.
-
|
@GiuliaIsNotAvailable just want to point out the changes were implemented a while ago, in xclim===0.48 , in case it wasn't clear. There was an important change in the API as of version 0.49 for SPI/SPEI, but the functions were improved as a result, I hope the documentation is clear enough for these changes, let me know if you have questions. Thanks again for helping us improve xclim. |
Beta Was this translation helpful? Give feedback.
Hi @aulemahal. I shall wait with confidence a code update. Thank you very much.