climate_indices
climate_indices copied to clipboard
monocongo
Can you help me understand how to fix the ValueError: Unable to calculate Pearson Type III parameters due to invalid L-moments
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
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.
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 https://github.com/monocongo/climate_indices/issues/512#issuecomment-1557373087, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJ5OUUYMEMIQXI6JBVTJS7LXHN5INANCNFSM6AAAAAAXVXSRYU . You are receiving this because you authored the thread.Message ID: @.***>
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.
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.
@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.
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
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
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.
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
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.
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 https://github.com/monocongo/climate_indices/issues/512#issuecomment-1615933063, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJ5OUU6RAFMIARDVT2YKSPDXOAVEVANCNFSM6AAAAAAXVXSRYU . You are receiving this because you authored the thread.Message ID: @.***>
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.
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.
@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 !
@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.
@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.
我在计算月尺度SPI Pearson 也碰到了同样的问题，可能是由于时间序列数据的问题，我将其加上很小的随机数就可计算。 nc_pre[yindex[i],xindex[i],:] = data + np.random.random(data.size)*0.01
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?
