Skip to content

SPEI : can't figure it out #1564

@GiuliaIsNotAvailable

Description

@GiuliaIsNotAvailable

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

Dataset.zip

Code of Conduct

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    supportQuestions and help for users/developers

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions