-
Notifications
You must be signed in to change notification settings - Fork 70
Description
Setup Information
- Xclim version: 0.46.0
- Python version: 3.10.13 (Jupyter notebook)
- Operating System: Windows 11
Context
Hello, 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:
import xclim.indices as indices
import xarray as xr
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']
spei_feb = indices.standardized_precipitation_evapotranspiration_index(wb = bic_daily_foramonth, freq = None, window = 1, dist = 'gamma', method = 'APP')
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:
agERA5_202302 = xr.open_dataset('202302_agERA5.nc').to_dataframe().reset_index()
agERA5_202302 = agERA5_202302.rename(columns={'time':'day'})
agERA5_202302 = agERA5_202302.set_index(['day','latitude','longitude'])
agERA5_202302['bic_daily'] = agERA5_202302['ppn'] - agERA5_202302['pev']
agERA5_202302 = agERA5_202302.to_xarray()
bic_daily_foramonth = agERA5_202302['bic_daily']
If I run:
spei_feb = atmos.standardized_precipitation_evapotranspiration_index(wb = bic_daily_foramonth, freq = None, window = 1, dist = 'gamma', method = 'APP')
I got: ValueError: Inputs have different frequencies. Got : [].To mute this, set xclim's option data_validation='log'.
If I run:
spei_feb = atmos.standardized_precipitation_evapotranspiration_index(wb = bic_daily_foramonth, freq = None, window = 1, dist = 'gamma', method = 'APP')
I got: KeyError: 'units'
I decided to try again with some other changes. Below the last python code:
agERA5_202302 = xr.open_dataset('202302_agERA5.nc').to_dataframe().reset_index()
agERA5_202302 = agERA5_202302.rename(columns={'time':'day'})
agERA5_202302 = agERA5_202302.set_index(['day','latitude','longitude'])
agERA5_202302['bic_daily'] = agERA5_202302['ppn'] - agERA5_202302['pev']
agERA5_202302 = agERA5_202302.to_xarray()
bic_daily_foramonth = agERA5_202302['bic_daily']
bic_daily_foramonth.attrs['units'] = 'mm / day'
spei_feb = indices.standardized_precipitation_evapotranspiration_index(wb = bic_daily_foramonth, freq = None, window = 1, dist = 'gamma', method = 'APP')
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
I'm sorry for being wordy. Any suggestions?
Thank you in advance
Giulia
Code of Conduct
- I agree to follow this project's Code of Conduct