Comparison with the R package

[GiuliaIsNotAvailable]

Dear @monocongo,
I'm part of a research group involved in agrometeorological studies and we have developed a service to calculate the SPEI index at a monthly frequency, applying the SPEI R package. Now we need to integrate this product with a daily estimate (with daily frequency), for which we would like to use your climate_indices Python package.
First of all, we encountered a problem because the SPEI R estimate is based on the log-logistic probability density function, whereas your package allows for choosing between the Pearson III and gamma probability density function. We wonder if you are planning to add some other functions, in particular the log-logistic one (as suggested by Vicente-Serrano) or if you noticed that this last one is worse than the others.
For the moment we have compared the SPEI values obtained from the R package with those produced by climate_indices, setting the Pearson III pdf for both computations (on the same dataset), as shown in the attached file SPEI_daily_Pearson.png. We noticed that R SPEI values are quite close to climate_indices values only during winter months while in the other months they are completely different and even reach up to '-inf '. Do you have any suggestions to understand this behaviour?
Below we show our R code for the calculation, where cwb stands for 'climatic water balance' (precipitation-evapotranspiration):

spei.cell <- spei(

ts(cwb, freq=365, start=c(1981,1,1)),

scale = 30,

kernel= list(type="rectangular", shift=0),

distribution="PearsonIII",

fit= "ub-pwm",

params=NULL,

na.rm=F,

ref.start=NULL,

ref.end=NULL

)

In attachment there is the Python code (SPEI_daily.txt). In addition we would like to know the parameter setting in your package to compare it with our setting in the R code (or how can we change them in order to set them as in the R code?). We leave also a sample of our data (ppn=precipitation, pev=evapotranspiration), in case you would try to test our code.

Thank you in advance!

E5L_8810_ppn_pev_19812021.csv

[monocongo]

Thanks for submitting this issue @GiuliaIsNotAvailable

I was not able to run your example code as-is, but I was able to get a result using the code below:

import pandas as pd

from climate_indices import compute, indices

data_csv = "/Users/jadams/tmp/E5L_8810_ppn_pev_19812021.csv"
df = pd.read_csv(filepath_or_buffer=data_csv, index_col=0)

df.index = pd.to_datetime(df.index)

prec_rsum = df["ppn"].rolling("30D", min_periods=30).sum().dropna()
evap_rsum = df["pev"].rolling("30D", min_periods=30).sum().dropna()

initial_year = prec_rsum.index[0].year
calibration_year_initial = prec_rsum.index[0].year
calibration_year_final = prec_rsum.index[-1].year
period_times = 366
scale = 1
periodicity = compute.Periodicity.daily

speival = indices.spei(
    precips_mm= prec_rsum.values,
    pet_mm= evap_rsum.values,
    scale=scale,
    distribution=indices.Distribution.pearson,
    data_start_year=initial_year,
    calibration_year_initial=calibration_year_initial,
    calibration_year_final=calibration_year_final,
    periodicity=compute.Periodicity.daily,
)

I ran this code using the latest dev branch (which will be the new main branch once I merge #526) and it produced reasonable-looking values. I have no idea as to what the R package is doing nor how to run it to get corresponding data.

You might want to collaborate with the team at @Ouranosinc/xclim which has put together an even better package for climate indices than this one, and they even have a current issue where they're developing some SPEI code that might be of interest to your group. This one may be of interest as well: Ouranosinc/xclim#1474

[monocongo]

@GiuliaIsNotAvailable You may want to contact the group who wrote the R code here: https://github.com/sbegueria/SPEI

On their page, they mention that the package fetched from the CRAN repo may not be the latest code, so you might get some mileage out of trying again with their latest R code from GitHub instead.

Suerte!

[GiuliaIsNotAvailable]

@monocongo thanks you for your advices and assistance!