Can you help me understand how to fix the ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments

Question

Can you help me understand how to fix the ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments

Opened this issue a year ago · 23 comments

I had a problem with L-moments when i excute my climate indices command:
rocess_climate_indices --index spi --periodicity monthly --netcdf_precip sliced_pre_1901_18.nc --var_name_precip pre --output_file_base sliced_pr.nc --scales 12 --calibration_start_year 1901 --calibration_end_year 2018
What's the problem ! @monocongo @monocongo @WeatherGod @ScottWales @dawiedotcom

Answer 1 · 2023-05-22T14:58:03.000Z

Please post a link to your datasets so I can attempt to reproduce and understand the error. For example sliced_pre_1901_18.nc There have been many issues regarding L-moments in the past, and all seem to have been resolved now if using the latest scipy etc. Please install the climate_indices package from the master branch (which is ahead of what's now on PyPI) to see if that will fix your issue. My apologies, I haven't had time to update PyPI with the latest/greatest code yet.

Answer 2 · 2023-06-05T07:54:53.000Z

I have now an output for spi but with masked array, when I visualize the nc file there are a warning in my data masked array (warning valid max and min). I'll send you my output spi to check if my output is correct or not. Are you recieve my output drive link ?

…

On Mon, May 22, 2023, 15:58 James Adams ***@***.***> wrote: Please post a link to your datasets so I can attempt to reproduce and understand the error. For example sliced_pre_1901_18.nc There have been many issues regarding L-moments in the past, and all seem to have been resolved now if using the latest scipy etc. Please install the climate_indices package from the master branch (which is ahead of what's now on PyPI) to see if that will fix your issue. My apologies, I haven't had time to update PyPI with the latest/greatest code yet. — Reply to this email directly, view it on GitHub <#512 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJ5OUUYMEMIQXI6JBVTJS7LXHN5INANCNFSM6AAAAAAXVXSRYU> . You are receiving this because you authored the thread.Message ID: ***@***.***>

Answer 3 · 2023-06-07T14:30:54.000Z

Please post a link to your data here.

Are you using the latest code from the master branch of this repository, or are you using the version from PyPI which is a little behind that (sorry for out of sync, tech debt)? Please use the latest code from this repository which should have the L-moments issue worked out. There have been many problems with the L-moments calculations over the years but seems like they've been rectified now.

Answer 4 · 2023-06-13T06:24:24.000Z

I'm having a similar issue, which I've managed to trace to a violation of theorem 1 in Hosking (1990), which bounds tau3 (or t3 in the code here) to -1 < t3 < 1. In my case, I get t3 = 1. Below is a summary of how I found it.

I have precipitation and evapotranspiration data with shape (lat, lon, time) = (8520, 7320, 180), in monthly time steps from Jan 2004 to Dec 2018. I'm trying to calculate SPEI.

I added some print messages to the main loop in main._apply_along_axis_double() to get the lat, lon indexes of where the L-moments error were happening. Then I proceeded with a step-by-step calculation using the modules. Below is an example of what is causing the problem.

For the month of August, we have:

precipitation = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
evapotranspiration = [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]
time_step_values = [2997., 2997., 2997., 2997., 2997., 2994., 2997., 2997., 2997., 2997., 2997., 2997., 2997., 2997., 2997.]
lmoments = [ 2.9968e+03  2.0000e-01 -1.0000e+00]

When the calculation calls lmoments._estimate_pearson3_parameters(), we get

t3 = abs(lmoments[2])  # L-skewness?
t3 >= 1
True

And hence, the error. Using scipy.stats.pearson3.fit() on the data with both maximum likelihood or method of moments (MM) we get valid output:

scipy.stats.pearson3.fit(time_step_values, method='MM')
(-3.474284756042035, 2996.800000000001, 0.7483314723274126)

Just for comparison, for the month of February, which doesn't trigger an error, we get estimates that are pretty close between scipy with "MM" and the estimate_pearson3_parameters routine:

scipy: {'loc': 3008.3333333333703, 'skew': 0.5728926162569001, 'scale': 4.422166348050119}
estimate_pearson3_parameters: {'loc': 3008.3333333333335, 'skew': 0.6312923338011005, 'scale': 4.68319994913074}

Looking at other points this seems to happen when all P = 0 for a certain month (where seasonality is high). Maximum likelihood has been used before to calculate SPI/SPEI, with L-moments providing initial values (as in Stagge et al. 2015, https://doi.org/10.1002/joc.4267), so this might be an option in the future for similar cases.

Answer 5 · 2023-06-13T16:24:37.000Z

@caiomattos' analysis looks really solid here. Thanks for the input! If there's a way to fix our code to accommodate this then a PR will be welcome.

Answer 6 · 2023-06-21T09:38:12.000Z

Hi,

I installed climate-indices into a fresh Python3.10 virtual env with the latest scipy version 1.10.1 and climate indices 1.10.13
and am getting the same error:
Traceback (most recent call last):
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 125, in worker
result = (True, func(*args, **kwds))
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
return list(map(*args))
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1428, in _apply_along_axis_double
computed_array[i, j] = func1d(x[j], y[j], parameters=params["args"])
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1151, in _spei
return indices.spei(precips_mm=precips,
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/indices.py", line 333, in spei
compute.transform_fitted_pearson(
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/compute.py", line 553, in transform_fitted_pearson
pearson_parameters(
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/compute.py", line 349, in pearson_parameters
params = lmoments.fit(time_step_values)
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/lmoments.py", line 30, in fit
raise ValueError(message)
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1782, in main
_compute_write_index(kwrgs)
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1044, in _compute_write_index
_parallel_process(keyword_arguments["index"],
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1324, in _parallel_process
pool.map(_apply_along_axis_double, chunk_params)
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 364, in map
return self._map_async(func, iterable, mapstar, chunksize).get()
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 771, in get
raise self._value
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
multiprocessing.pool.RemoteTraceback:
"""
Traceback (most recent call last):
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 125, in worker
result = (True, func(*args, **kwds))
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
return list(map(*args))
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1428, in _apply_along_axis_double
computed_array[i, j] = func1d(x[j], y[j], parameters=params["args"])
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1151, in _spei
return indices.spei(precips_mm=precips,
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/indices.py", line 333, in spei
compute.transform_fitted_pearson(
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/compute.py", line 553, in transform_fitted_pearson
pearson_parameters(
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/compute.py", line 349, in pearson_parameters
params = lmoments.fit(time_step_values)
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/lmoments.py", line 30, in fit
raise ValueError(message)
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
File "/home/dray/anaconda3/envs/CI/bin/process_climate_indices", line 8, in
sys.exit(main())
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1782, in main
_compute_write_index(kwrgs)
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1044, in _compute_write_index
_parallel_process(keyword_arguments["index"],
File "/home/dray/anaconda3/envs/CI/lib/python3.10/site-packages/climate_indices/main.py", line 1324, in _parallel_process
pool.map(_apply_along_axis_double, chunk_params)
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 364, in map
return self._map_async(func, iterable, mapstar, chunksize).get()
File "/home/dray/anaconda3/envs/CI/lib/python3.10/multiprocessing/pool.py", line 771, in get
raise self._value
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments

On running:
process_climate_indices --index spei --periodicity monthly --netcdf_precip LME_monthlyRAIN_AUSTNZ.nc --var_name_precip RAIN --netcdf_pet LME_monthlyEto_AUSTNZ.nc --var_name_pet PM_FAO_56 --output_file_base nclimgrid_SPEI12 --scales 12 --calibration_start_year 1986 --calibration_end_year 2005 --multiprocessing all

Link to rainfall and PET datasets:
https://www.dropbox.com/scl/fi/8vpxhlyb7bllihlfm0th7/LME_monthlyRAIN_AUSTNZ.nc?dl=0&rlkey=wa2egqkal8oxpcmugyb46b647
https://www.dropbox.com/scl/fi/3kot9u44qpdt5bgp7v14p/LME_monthlyEto_AUSTNZ.nc?dl=0&rlkey=4uj8pgywhs1ffckuo3rrhqp21

Not sure if related or not but also get:
NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.

I have to run : export NUMBA_DISABLE_JIT=1 to get the job to even start to process to get to the L-moments error

I got the above command to run and worked a couple of weeks ago using version 1.10.12 but it wouldnt run this time so I updated to 1.10.13 and still not working.

conda list:
_libgcc_mutex 0.1 conda_forge conda-forge
_openmp_mutex 4.5 2_gnu conda-forge
blosc 1.21.4 h0f2a231_0 conda-forge
bzip2 1.0.8 h7f98852_4 conda-forge
c-ares 1.19.1 hd590300_0 conda-forge
ca-certificates 2023.5.7 hbcca054_0 conda-forge
certifi 2023.5.7 pyhd8ed1ab_0 conda-forge
cftime 1.6.2 py310hde88566_1 conda-forge
click 8.1.3 pypi_0 pypi
climate-indices 1.0.13 pypi_0 pypi
cloudpickle 2.2.1 pypi_0 pypi
curl 8.1.2 h409715c_0 conda-forge
dask 2023.6.0 pypi_0 pypi
esmf 8.4.2 nompi_h20110ff_0 conda-forge
expat 2.5.0 hcb278e6_1 conda-forge
fsspec 2023.6.0 pypi_0 pypi
gsl 2.7 he838d99_0 conda-forge
hdf4 4.2.15 h501b40f_6 conda-forge
hdf5 1.14.0 nompi_hb72d44e_103 conda-forge
icu 72.1 hcb278e6_0 conda-forge
importlib-metadata 6.7.0 pypi_0 pypi
keyutils 1.6.1 h166bdaf_0 conda-forge
krb5 1.20.1 h81ceb04_0 conda-forge
ld_impl_linux-64 2.40 h41732ed_0 conda-forge
libaec 1.0.6 hcb278e6_1 conda-forge
libblas 3.9.0 17_linux64_openblas conda-forge
libcblas 3.9.0 17_linux64_openblas conda-forge
libcurl 8.1.2 h409715c_0 conda-forge
libedit 3.1.20191231 he28a2e2_2 conda-forge
libev 4.33 h516909a_1 conda-forge
libexpat 2.5.0 hcb278e6_1 conda-forge
libffi 3.4.2 h7f98852_5 conda-forge
libgcc-ng 13.1.0 he5830b7_0 conda-forge
libgfortran-ng 13.1.0 h69a702a_0 conda-forge
libgfortran5 13.1.0 h15d22d2_0 conda-forge
libgomp 13.1.0 he5830b7_0 conda-forge
libiconv 1.17 h166bdaf_0 conda-forge
libjpeg-turbo 2.1.5.1 h0b41bf4_0 conda-forge
liblapack 3.9.0 17_linux64_openblas conda-forge
libnetcdf 4.9.2 nompi_h0f3d0bb_105 conda-forge
libnghttp2 1.52.0 h61bc06f_0 conda-forge
libnsl 2.0.0 h7f98852_0 conda-forge
libopenblas 0.3.23 pthreads_h80387f5_0 conda-forge
libsqlite 3.42.0 h2797004_0 conda-forge
libssh2 1.11.0 h0841786_0 conda-forge
libstdcxx-ng 13.1.0 hfd8a6a1_0 conda-forge
libuuid 2.38.1 h0b41bf4_0 conda-forge
libxml2 2.11.4 h0d562d8_0 conda-forge
libzip 1.9.2 hc929e4a_1 conda-forge
libzlib 1.2.13 hd590300_5 conda-forge
llvmlite 0.40.1rc1 pypi_0 pypi
locket 1.0.0 pypi_0 pypi
lz4-c 1.9.4 hcb278e6_0 conda-forge
nco 5.1.6 hd62b316_0 conda-forge
ncurses 6.4 hcb278e6_0 conda-forge
netcdf-fortran 4.6.1 nompi_h4f3791c_100 conda-forge
netcdf4 1.6.4 nompi_py310hde23a83_100 conda-forge
numba 0.57.0 pypi_0 pypi
numpy 1.24.3 pypi_0 pypi
openssl 3.1.1 hd590300_1 conda-forge
packaging 23.1 pypi_0 pypi
pandas 2.0.2 pypi_0 pypi
partd 1.4.0 pypi_0 pypi
pip 23.1.2 pyhd8ed1ab_0 conda-forge
python 3.10.0 h543edf9_3_cpython conda-forge
python-dateutil 2.8.2 pypi_0 pypi
python_abi 3.10 3_cp310 conda-forge
pytz 2023.3 pypi_0 pypi
pyyaml 6.0 pypi_0 pypi
readline 8.2 h8228510_1 conda-forge
scipy 1.10.1 pypi_0 pypi
setuptools 67.7.2 pyhd8ed1ab_0 conda-forge
six 1.16.0 pypi_0 pypi
snappy 1.1.10 h9fff704_0 conda-forge
sqlite 3.42.0 h2c6b66d_0 conda-forge
tempest-remap 2.2.0 h43474b4_0 conda-forge
tk 8.6.12 h27826a3_0 conda-forge
toolz 0.12.0 pypi_0 pypi
tzdata 2023.3 pypi_0 pypi
udunits2 2.2.28 hc3e0081_0 conda-forge
wheel 0.40.0 pyhd8ed1ab_0 conda-forge
xarray 2023.5.0 pypi_0 pypi
xz 5.2.6 h166bdaf_0 conda-forge
zipp 3.15.0 pypi_0 pypi
zstd 1.5.2 h3eb15da_6 conda-forge

Answer 7 · 2023-06-30T18:43:35.000Z

There is an imminent update with code that I've tested against these Australian datasets, looking good now using the development branch named issue_522_pyproject_poetry in case you want to try it out right away.

Answer 8 · 2023-07-01T06:31:31.000Z

Hi...thanks for that.

I did an uninstall and reinstall from the issue_522_pyproject_poetry tree.

The line 10 import typing in main.py needs to import Dict and Any as well.

When I fixed that, it ran, but again threw the L-moments error on Pearson typeIII :(

Is there a way to just run the Gamma distribution by the way?

And, from your answer, I assume the package is okay with cftime conventions?
I have run the example .nc grids and they work fine by the way

Answer 9 · 2023-07-01T14:05:51.000Z

Thanks @Darren-Ray I appreciate the beta testing!

On my copy of that branch, I am having good results running the processing command. I've created a fresh conda environment and installed from source, like so:

conda create -n myvenv poetry pytest
conda activate myvenv
python -m poetry install
python -m poetry run pytest
cd ~/tmp/data
process_climate_indices --index spei --periodicity monthly --netcdf_precip LME_monthlyRAIN_AUSTNZ.nc --var_name_precip RAIN --netcdf_pet LME_monthlyEto_AUSTNZ.nc --var_name_pet PM_FAO_56 --output_file_base nclimgrid_SPEI12 --scales 12 --calibration_start_year 1986 --calibration_end_year 2005 --multiprocessing all

I have added cftime as a dependency which got me past an error I saw around that. So we should be OK for cftime conventions, but I'm not well-versed in that and can't promise too much there.

We used to have an option that allowed users to specify a single distribution but looking at the docs it seems that that was removed at some point, or maybe I am misremembering. Anyway, it shouldn't be too hard to re-implement that by adding a --distfit option or something like that.

Answer 10 · 2023-07-10T15:10:06.000Z

Hello every one, I tested the package "Another one" with fixing some problems and it's working fine. I would thank all of you, specially Mr. James Adams. All my good wishes to you. Le sam. 1 juil. 2023 à 15:06, James Adams ***@***.***> a écrit :

…

Thanks @Darren-Ray <https://github.com/Darren-Ray> I appreciate the beta testing! On my copy of that branch, I am having good results running the processing command. I've created a fresh conda environment and installed from source, like so: conda create -n myvenv poetry pytest conda activate myvenv python -m poetry install python -m poetry run pytest cd ~/tmp/data process_climate_indices --index spei --periodicity monthly --netcdf_precip LME_monthlyRAIN_AUSTNZ.nc --var_name_precip RAIN --netcdf_pet LME_monthlyEto_AUSTNZ.nc --var_name_pet PM_FAO_56 --output_file_base nclimgrid_SPEI12 --scales 12 --calibration_start_year 1986 --calibration_end_year 2005 --multiprocessing all I have added cftime as a dependency which got me past an error I saw around that. So we should be OK for cftime conventions, but I'm not well-versed in that and can't promise too much there. We used to have an option that allowed users to specify a single distribution but looking at the docs it seems that that was removed at some point, or maybe I am misremembering. Anyway, it shouldn't be too hard to re-implement that by adding a --distfit option or something like that. — Reply to this email directly, view it on GitHub <#512 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJ5OUU6RAFMIARDVT2YKSPDXOAVEVANCNFSM6AAAAAAXVXSRYU> . You are receiving this because you authored the thread.Message ID: ***@***.***>

Answer 11 · 2024-01-14T05:08:53.000Z

Hi! @monocongo
I have the same error, ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments, when calculating spi using data from CRU (https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.07/cruts.2304141047.v4.07/pre/)
I try to install it from the master branch and update Scipy but it is still not working.
Any suggestion?

The example data from nclimgrid and cmorph is working well.
The data from CHIRPS and ERA5 is working after reordering the variables and editing attributes.

Answer 12 · 2024-01-22T18:56:17.000Z

Thanks, @chaowats

Please confirm that you get the same result if installing this package into a fresh virtual environment using poetry. Please install from PyPI rather than github.

Answer 13 · 2024-01-23T13:30:29.000Z

@chaowats You can tell me how you installed the climate indices package because I'm traveling with the CHIRPS data.
I installed the package and reorder the lat lon time, but it still not working !

Answer 14 · 2024-01-27T06:58:30.000Z

@monocongo Yes I did install to a new environment with the 2.0 version and master branch using pip install climate-indices from PyPI (but with conda environment, is that okay?)
I got both the same Pearson Type III error.

Answer 15 · 2024-01-27T07:02:37.000Z

@ben29med I install using pip install climate-indices from PyPI but with the CHIRPS data I change precipitation units from mm/month to mm, longitude to lon, latitude to lat, and reorder to lat,lon,time before process_climate_indices --index spi.

Answer 16 · 2024-02-11T10:37:37.000Z

@chaowats please, send me the nclimgrid dataset example. Thank's !

Answer 17 · 2024-04-03T03:56:11.000Z

我在计算月尺度SPI Pearson 也碰到了同样的问题，可能是由于时间序列数据的问题，我将其加上很小的随机数就可计算。
nc_pre[yindex[i],xindex[i],:] = data + np.random.random(data.size)*0.01

Answer 18 · 2024-04-05T23:34:46.000Z

Hi! @monocongo I have the same error, ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments, when calculating spi using data from CRU (https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.07/cruts.2304141047.v4.07/pre/) I try to install it from the master branch and update Scipy but it is still not working. Any suggestion?

The example data from nclimgrid and cmorph is working well. The data from CHIRPS and ERA5 is working after reordering the variables and editing attributes.

Hello, I came across the same problem with CRU data but ERA5 works fine. Did you find the solution for CRU data?

Answer 19 · 2024-08-29T19:48:33.000Z

I also experience the L-moments error. I use the UCSB CHIRPS data. This follows a number of runtime warnings "divide by zero encountered in divide" and " invalid value encountered in multiply", which may be related to the L-moments error? Full log below.
@monocongo is there a reference data set which should work?

The data can be found here: https://tubcloud.tu-berlin.de/s/4YcpA9RXLBXzka6

I am running the following commands:

 `$ conda create -n ci310 python=3.10
  $ conda activate ci310
  $ git clone https://github.com/monocongo/climate_indices.git
  $ cd climate_indices/
  $ poetry install
  $ process_climate_indices --index spi --periodicity monthly --netcdf_precip africa_chirps_1months_1981_2023.nc --var_name_precip precip --output_file_base CHIRPS --scales 3 --calibration_start_year 1981 --calibration_end_year 2023

Am I doing this right?
Many thanks for you help, much appreciated.
Below is the full terminal output.

2024-08-29  19:16:29 INFO Start time:    2024-08-29 19:16:29.147045
 2024-08-29  19:16:38 INFO Computing 3-month SPI/Pearson
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:455: RuntimeWarning: divide by zero encountered in divide
  alpha = 4.0 / (skew * skew)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:459: RuntimeWarning: invalid value encountered in multiply
  return loc - ((alpha * scale * skew) / 2.0)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:455: RuntimeWarning: divide by zero encountered in divide
  alpha = 4.0 / (skew * skew)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:459: RuntimeWarning: invalid value encountered in multiply
  return loc - ((alpha * scale * skew) / 2.0)
subfolder/miniconda3/envs/ci310/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py:2176: RuntimeWarning: invalid value encountered in divide
  x = np.asarray((x - loc)/scale, dtype=dtyp)
2024-08-29  19:21:48 ERROR Unable to calculate Pearson Type III parameters due to invalid L-moments
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:455: RuntimeWarning: divide by zero encountered in divide
  alpha = 4.0 / (skew * skew)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:459: RuntimeWarning: invalid value encountered in multiply
  return loc - ((alpha * scale * skew) / 2.0)
subfolder/miniconda3/envs/ci310/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py:2176: RuntimeWarning: invalid value encountered in divide
  x = np.asarray((x - loc)/scale, dtype=dtyp)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:455: RuntimeWarning: divide by zero encountered in divide
  alpha = 4.0 / (skew * skew)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:459: RuntimeWarning: invalid value encountered in multiply
  return loc - ((alpha * scale * skew) / 2.0)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:455: RuntimeWarning: divide by zero encountered in divide
  alpha = 4.0 / (skew * skew)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:459: RuntimeWarning: invalid value encountered in multiply
  return loc - ((alpha * scale * skew) / 2.0)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:455: RuntimeWarning: divide by zero encountered in divide
  alpha = 4.0 / (skew * skew)
subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py:459: RuntimeWarning: invalid value encountered in multiply
  return loc - ((alpha * scale * skew) / 2.0)
subfolder/miniconda3/envs/ci310/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py:2176: RuntimeWarning: invalid value encountered in divide
  x = np.asarray((x - loc)/scale, dtype=dtyp)
2024-08-29  19:31:17 ERROR Unable to calculate Pearson Type III parameters due to invalid L-moments
2024-08-29  19:31:17 ERROR Failed to complete
multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
    return list(map(*args))
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1292, in _apply_along_axis
    computed_array = np.apply_along_axis(func1d, axis=axis_index, arr=sub_array, parameters=args)
  File "<__array_function__ internals>", line 200, in apply_along_axis
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/site-packages/numpy/lib/shape_base.py", line 402, in apply_along_axis
    buff[ind] = asanyarray(func1d(inarr_view[ind], *args, **kwargs))
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1102, in _spi
    return indices.spi(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/indices.py", line 184, in spi
    values = compute.transform_fitted_pearson(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py", line 624, in transform_fitted_pearson
    probabilities_of_zero, locs, scales, skews = pearson_parameters(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py", line 307, in pearson_parameters
    prob, loc, scale, skew = calculate_time_step_params(time_step_values)
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py", line 256, in calculate_time_step_params
    params = lmoments.fit(time_step_values)
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/lmoments.py", line 33, in fit
    raise ValueError(message)
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1617, in main
    _compute_write_index(kwrgs)
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 997, in _compute_write_index
    _parallel_process(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1256, in _parallel_process
    pool.map(_apply_along_axis, chunk_params)
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 367, in map
    return self._map_async(func, iterable, mapstar, chunksize).get()
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 774, in get
    raise self._value
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
    return list(map(*args))
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1292, in _apply_along_axis
    computed_array = np.apply_along_axis(func1d, axis=axis_index, arr=sub_array, parameters=args)
  File "<__array_function__ internals>", line 200, in apply_along_axis
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/site-packages/numpy/lib/shape_base.py", line 402, in apply_along_axis
    buff[ind] = asanyarray(func1d(inarr_view[ind], *args, **kwargs))
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1102, in _spi
    return indices.spi(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/indices.py", line 184, in spi
    values = compute.transform_fitted_pearson(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py", line 624, in transform_fitted_pearson
    probabilities_of_zero, locs, scales, skews = pearson_parameters(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py", line 307, in pearson_parameters
    prob, loc, scale, skew = calculate_time_step_params(time_step_values)
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/compute.py", line 256, in calculate_time_step_params
    params = lmoments.fit(time_step_values)
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/lmoments.py", line 33, in fit
    raise ValueError(message)
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "subfolder/miniconda3/envs/ci310/bin/process_climate_indices", line 6, in <module>
    sys.exit(main())
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1617, in main
    _compute_write_index(kwrgs)
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 997, in _compute_write_index
    _parallel_process(
  File "subfolder/SPI/climate_indices/climate_indices/climate_indices/climate_indices/src/climate_indices/__main__.py", line 1256, in _parallel_process
    pool.map(_apply_along_axis, chunk_params)
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 367, in map
    return self._map_async(func, iterable, mapstar, chunksize).get()
  File "subfolder/miniconda3/envs/ci310/lib/python3.10/multiprocessing/pool.py", line 774, in get
    raise self._value
ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments`

Answer 20 · 2024-09-01T15:14:22.000Z

@rabons I see the same error when I run on my laptop. Looking at the data with ncdump -h shows that the time values are using a periodicity of "dekad":

        float precip(lat, lon, time) ;
                precip:_FillValue = -9999.f ;
                precip:units = "mm" ;
                precip:standard_name = "convective precipitation rate" ;
                precip:long_name = "Climate Hazards group InfraRed Precipitation with Stations" ;
                precip:time_step = "dekad" ;
                precip:missing_value = -9999.f ;
                precip:geospatial_lat_min = "NaN" ;
                precip:geospatial_lat_max = "NaN" ;
                precip:geospatial_lon_min = "NaN" ;
                precip:geospatial_lon_max = "NaN" ;
                precip:grid_mapping = "crs" ;

I was only recently made aware of this time step period but apparently, it's used by some researchers, so I'd like to support it as a new Periodicity enumerated type. I suspect the issue concerns a mismatch between monthly and dekad time steps for this error case. For example, we have two period (time step) types we support -- monthly and daily, and this will require adding a third type and modifying the processing loops to account for 36.5 dekads per year. We needed to do this before to add support for daily data (various modifications to allow for 366-day years) so this will be another iteration of that kind of refactoring.

Answer 21 · 2024-09-01T18:23:52.000Z

@monocongo thank you for lookig into this!

time_step = "dekad" is actually a labeling error in this specific file, I think. Will double-check tomorrow. Would this label affect the SPI calculations in any way? I can re-produce it with a correct label.

EDIT: I have changed the file with correct labeling, it doesnt change the error with climate-indices. I have created this file with this script (slightly modified): https://bennyistanto.github.io/spi/chirpstif/

I actually do need a "dekadal" time step, but more important is that I make the monthly calculations work.

Answer 22 · 2024-09-02T16:18:57.000Z

OK then this isn't caused by the dekad time step since it's just a label. Thanks for the clarification @rabons

This error is raised in the lmoments.py code, here:

    # ensure the validity of the L-moments
    if (lmoments[1] <= 0) or (t3 >= 1):
        message = "Unable to calculate Pearson Type III parameters due to invalid L-moments"
        _logger.error(message)
        raise ValueError(message)

We can add some debugging around that to see if there are anomalies in the data etc. causing this to happen.

Answer 23 · 2024-09-05T09:35:04.000Z

Thank you.

EDIT: SPI calculations seem to be working for a file clipped to a shapefile (in this case, Uganda boundaries). There may have been an unrecognized problem with my previous input file. The working input file: https://tubcloud.tu-berlin.de/s/omStM55dAmC2XAw
and the spi result: https://tubcloud.tu-berlin.de/s/BE9E3ErefZD2MTn

The issue is hence solved for me when first clipping the input file (CHIRPS africa_monthly tiffs). Not sure what this is about. A few 'division by zero' errors persist, but I can still produce valid SPI output.