ausclimateservice / indices Goto Github PK

Code for calculating climate indices

Python 100.00%

indices's Introduction

indices

This directory contains two command line programs for calculating climate indices.

run_icclim.py uses the icclim library to calculate simple indices such as those used by the ETCCDI and Climdex projects. Detailed descriptions of each index are available in the icclim API and an associated ATBD.

run_xclim.py uses the xclim library to calculate a series of climate indicators documented at the xclim climate indicators API.

The major difference between the icclim and xclim indices is that the icclim indices represent a set of widely used, rigidly defined metrics (e.g. the RX5day index can't be modified to calculate the maximum rainfall over 4 days instead) while the xclim indicators can be configured to set different thresholds, times of year, etc.

Python environment

If you're a member of the xv83 project on NCI (i.e. if you're part of the Australian Climate Service), the easiest way to execute run_icclim.py is to run it at the command line with the version of Python installed at: /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python. For example:

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_icclim.py -h

If you'd like to run the script/s in your own Python environment, you'll need to install the following libraries using conda:

$ conda install -c conda-forge icclim cmdline_provenance gitpython

Usage: `run_icclim.py`

Running the script at the command line with the -h option explains the user options, including the very long list of indices that the script can calculate.

$ python run_icclim.py -h

Example 1: Simple indices

The most basic use of run_icclim.py requires passing the script the name of the index to calculate, the name of the output file, the name/s of the input file/s, and the name of the variable to access from that file/s, For example:

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_icclim.py TXx /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/AUS-15/BOM/ECMWF-ERA5/evaluation/r1i1p1f1/BOM-BARPA-R/v1/climdex/txx/txx_AUS-15_ECMWF-ERA5_evaluation_r1i1p1f1_BOM-BARPA-R_v1_year_197901-200112.nc --input_files /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/output/AUS-15/BOM/ECMWF-ERA5/evaluation/r1i1p1f1/BOM-BARPA-R/v1/day/tasmax/tasmax_AUS-15_ECMWF-ERA5_evaluation_r1i1p1f1_BOM-BARPA-R_v1_day_*.nc --variable tasmax --verbose

In the example above, the --verbose flag has also been invoked so that the program prints its progress to the screen.

Example 2: Sub-daily input data

All the icclim climate indices are calculated from daily data. If you want to input sub-daily (e.g. hourly) data instead, you need to use the --time_agg option to specify how to temporally aggregate the data. For instance,

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_icclim.py TNn /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/GLOBAL-gn/none/ECMWF-ERA5/evaluation/r1i1p1f1/none/none/climdex/tnn/tnn_AUS-gn_ECMWF-ERA5_evaluation_r1i1p1f1_year_195901-202112.nc --input_files /g/data/rt52/era5/single-levels/reanalysis/2t/*/*.nc --variable t2m --time_agg min --start_date 1959-01-01 --end_date 2021-12-31 --lon_bnds 111.975 156.275 --lat_bnds -44.525 -9.975 --hshift --verbose

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_icclim.py R10mm /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/AUS-gn/none/ECMWF-ERA5/evaluation/r1i1p1f1/none/none/climdex/r10mm/r10mm_AUS-gn_ECMWF-ERA5_evaluation_r1i1p1f1_year_195901-202112.nc --input_files /g/data/rt52/era5/single-levels/reanalysis/mtpr/*/*.nc --variable mtpr --time_agg mean --start_date 1959-01-01 --end_date 2021-12-31 --lon_bnds 111.975 156.275 --lat_bnds -44.525 -9.975 --hshift --verbose

In both these examples, the --hshift option has been used to shift the time axis values back an hour (which is required when working with ERA5 data). The --start_date, --end_date, --lat_bnds and --lon_bnds options have also been used to subset the data. Selecting a latitude / longitude box around Australia rather than processing the whole globe (this example uses the AGCD dataset spatial bounds) can be particularly useful. If the TNn example above is run without the --lat_bnds and --lon_bnds selection it takes several hours to process the whole ERA5 global grid (when given 120GB of RAM on the Gadi "hugemem" job queue). With the lat/lon subsetting it takes less than an hour.

Example 3: Multi-variate indices

For a bivariate index like dtr (mean diurnal temperature range), you need to use invoke the --input_files and --variable options twice. The first time to specify the maximum temperature files and variable name, and the second time to specify the minimum temperature files and variable name. For example:

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_icclim.py DTR /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/AUS-r005/none/BOM-AGCD/historical/v1/none/none/climdex/dtr/dtr_AUS-r005_BOM-AGCD_historical_v1_year_191001-202112.nc  --input_files /g/data/xv83/agcd-csiro/tmax/daily/tmax_AGCD-CSIRO_r005_*.nc --variable tmax --input_files /g/data/xv83/agcd-csiro/tmin/daily/tmin_AGCD-CSIRO_r005_*.nc --variable tmin --start_date 1910-01-01 --end_date 2021-12-31 --verbose

CF-compliance

The icclim library (because of the xclim library it builds upon) requires CF-compliant variable attributes (see the health checks and metadata attributes section of the xclim docs for details). The run_icclim.py script has metadata fixing capability built in that should work for most datasets so that users don't have to edit their files before passing them to run_icclim.py. If you pass the script a file that does cause a metadata-related error, you can start a new issue here describing the error and we can add appropriate metadata handling to fix the problem.

Example 4: "Along the time axis" indices

A small number of climate indices require calculations to be performed along the entire time axis. For instance, in the example below the R95pTOT index requires calculation of the 95th percentile along the time axis for the entire 1900-2021 period:

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_icclim.py R95pTOT /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/AUS-r005/none/BOM-AGCD/historical/v1/none/none/climdex/r95ptot/r95ptot_AUS-r005_BOM-AGCD_historical_v1_year_190001-202112.nc --input_files /g/data/xv83/agcd-csiro/precip/daily/precip-total_AGCD-CSIRO_r005_19000101-20220405_daily_space-chunked.zarr --variable precip --start_date 1900-01-01 --end_date 2021-12-31 --verbose

Indices that involve calculations along the entire time axis are much more memory intensive. The example above used about 115 GB of RAM on the Australian Research Environment at NCI and took about 15 minutes to run.

The icclim documentation has extensive notes on how to improve performance for indices like this via data chunking and parallelisation using dask. The run_icclim.py script has a --local_cluster option that can be used (in conjunction with the --nworkers and --nthreads options) to launch and configure a local dask cluster. We are still figuring out the optimal local cluster settings for different use cases (and will add advice to this documentation if we find any useful approaches), but in general the computation will be fastest if you simply allocate a large amount of memory to the job (e.g. using the largemem job queue on NCI) as opposed to fiddling around with a local dask cluster. The only time allocating a large amount of memory possibly won't be sufficient is when using the --base_period option, which seems to dramitically increase compute requirements. We haven't been able to find dask cluster settings to overcome this problem, so a --nslices option is available to slice up the input data (along the longitude axis) and process each slice one-by-one before putting it all back together at the end.

Usage: `run_xclim.py`

Running the script at the command line with the -h option explains the user options, including the list of indices that the script can calculate.

$ python run_xclim.py -h

Example 1: Simple indices

The most basic use of run_xclim.py requires passing the script the name of the index to calculate, the name of the output file, the name/s of the input file/s, the name of the variable to access from that file/s, a threshold value for the index, and date bounds if analysing a restricted part of the year. For example:

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_xclim.py frost_days frost_days.nc --input_files /g/data/wp00/data/QQ-CMIP5/ACCESS1-0/tasmin/rcp45/2036-2065/tasmin_AUS_ACCESS1-0_rcp45_r1i1p1_CSIRO-QQS-AGCD-1981-2010_day_wrt_1986-2005_2036-2065.nc --variable tasmin --thresh "0.1 degC" --date_bounds 08-01 10-15 --verbose

In the example above, the --verbose flag has also been invoked so that the program prints its progress to the screen.

Example 2: Indices involving daily mean temperature

Many datasets archive daily minimum and maximum temperature, but not daily mean temperature. For metrics that require daily mean temperature, run_xclim.py will calculate the mean if you input tasmin and tasmax. e.g.

$ /g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python run_xclim.py growing_degree_days gdd.nc --input_files /g/data/wp00/data/QQ-CMIP5/ACCESS1-0/tasmax/rcp45/2036-2065/tasmax_AUS_ACCESS1-0_rcp45_r1i1p1_CSIRO-QQS-AGCD-1981-2010_day_wrt_1986-2005_2036-2065.nc --variable tasmax --input_files /g/data/wp00/data/QQ-CMIP5/ACCESS1-0/tasmin/rcp45/2036-2065/tasmin_AUS_ACCESS1-0_rcp45_r1i1p1_CSIRO-QQS-AGCD-1981-2010_day_wrt_1986-2005_2036-2065.nc --variable tasmin --thresh "0 degC" --date_bounds 04-01 10-31 --verbose

The trick here was to use the --input_files and --variable options twice; once for the tasmax data and once for the tasmin.

Common configurations

Late season heat risk / wheat heat risk

Description: Number of days per year where the maximum temperature is > 32C between 1 Aug and 30 Nov
Arguments: tx_days_above metric with --thresh "32 degC" and --date_bounds 08-01 11-30

Heat risk fertility metric / sheep heat risk

Description: Number of days per year where the maximum temperature is > 32C between 15 Jan and 15 Jun
Arguments: tx_days_above metric with --thresh "32 degC" and --date_bounds 01-15 06-15

indices's People

Contributors

Stargazers

Watchers

Forkers

wendyksharples wcarthur andrew-dowdy ngben mitchellblack akumenyi

indices's Issues

TX90p/TX10p/TN90p/TN10p calculations do not complete when --base_period argument is used

Using the --base_period argument for percentiles seems to hang with AGCD tasmax/tasmin data and the job runs out of walltime. I haven't applied it yet to other models. Below is the output from one of failed jobs.

/g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python /g/data/xv83/users/bxn599/ACS/icclim/run_icclim.py --slice_mode year --verbose --start_date 1979-01-01 --end_date 2021-12-31 --base_period 1985-01-01 2014-12-31 --input_files [AGCD tasmax files on xv83] --variable tmax --drop_time_bounds TX10p TX10p_AUS-r005_BOM-AGCD_historical_none_year_19790101-20211231.nc
[########################################] | 100% Completed | 10m 17s

rm: cannot remove 'TX10p_AUS-r005_BOM-AGCD_historical_none_year_19790101-20211231.nc': No such file or directory
2023-03-11 09:06:12,661 Array size: (15706, 691, 886)
2023-03-11 09:06:12,662 Chunk size: Frozen({'time': (15706,), 'lat': (35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 26), 'lon': (60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 46)})
2023-03-11 09:06:12,662    ********************************************************************************************
2023-03-11 09:06:12,662    *                                                                                          *
2023-03-11 09:06:12,662    *          icclim                6.1.3   *
2023-03-11 09:06:12,662    *                                                                                          *
2023-03-11 09:06:12,662    *                                                                                          *
2023-03-11 09:06:12,662    *          Fri Mar 10 22:06:12 2023                                                    *
2023-03-11 09:06:12,662    *                                                                                          *
2023-03-11 09:06:12,662    *          BEGIN EXECUTION                                                                 *
2023-03-11 09:06:12,662    *                                                                                          *
2023-03-11 09:06:12,662    ********************************************************************************************
2023-03-11 09:06:12,662 Processing: 0%
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/dataset.py:4880: PerformanceWarning: Reshaping is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array.reshape(shape)

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array.reshape(shape)Explictly passing ``limit`` to ``reshape`` will also silence this warning
    >>> array.reshape(shape, limit='128 MiB')
  result = result._unstack_full_reindex(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/variable.py:1722: PerformanceWarning: Reshaping is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array.reshape(shape)

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array.reshape(shape)Explictly passing ``limit`` to ``reshape`` will also silence this warning
    >>> array.reshape(shape, limit='128 MiB')
  result = result._stack_once(dims, new_dim)
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xclim/core/cfchecks.py:44: UserWarning: Variable does not have a `cell_methods` attribute.
  _check_cell_methods(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xclim/core/cfchecks.py:48: UserWarning: Variable does not have a `standard_name` attribute.
  check_valid(vardata, "standard_name", data["standard_name"])
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 15
  result = blockwise(
=>> PBS: job killed: walltime 172826 exceeded limit 172800

Add support for save_thresholds parameter

Hi Damien, if it's not too much trouble can you add support for the save_thresholds parameter?

This option saves the threshold values to the output netCDF file when calculating percentiles. It will be useful to compare historical/future thresholds and also for calculating future percentiles relative to the historical threshold. Thanks!

Temporary fix to calculate tas indices for BARPA

Hi Damien,

Sorry I don't know how to create a pull request in git. Could you please modify L220 in run_icclim.py to:
elif variable in ['tas', 't2m', '2t', 'tasmean']:
Including tasmean will support BARPA runs which have currently named tas as tasmean.

Thanks!

Custom indices

Custom indices could be added to the script using:
https://icclim.readthedocs.io/en/stable/references/custom_indices.html#custom-indices

lon bnds selection on ERA5 data is not working properly

I am trying to calculate indices in ERA5 using rt52 data and the whole BARPA domain (54S-14N, 88E-208E). If I put these lat lon bounds into run_icclim.py the output domain is not correct. I also tried lon_bnds 88.0 -152.0 which gives the same result. Could this issue be caused by ERA5 being split at 180 and spanning -180/180?

ncdump -c era5_rt52/GLO-r025/none/ERA5/historical/none/none/v1/climdex/tx/tx_GLO-r025_ERA5_historical_v1_year_19950101-20141231.nc output gives the following (cut for clarity):

/g/data/xv83/bxn599/ACS/data/era5/raw_test/mx2t/mx2t_era5_oper_sfc_20220701-20220731.nc --variable mx2t --lat_bnds -54.0 14.0 --lon_bnds 88.0 208.0 --time_agg max --hshift --drop_time_bounds tx /g/data/xv83/bxn599/ACS/icclim_indices/era5_rt52/GLO-r025/none/ERA5/historical/none/none/v1/climdex/tx/tx_GLO-r025_ERA5_historical_v1_year_19950101-20141231.nc"

lon = -180, -179.75, -179.5, -179.25, -179, -178.75, -178.5, -178.25, -178, -177.75, -177.5, -177.25, -177, -176.75, -176.5, -176.25, -176, -175.75, -175.5, -175.25, -175, -174.75, -174.5, -174.25, -174, -173.75, -173.5, -173.25, -173, -172.75, -172.5, -172.25, -172, -171.75, -171.5, -171.25, -171, -170.75, -170.5, -170.25, -170, -169.75, -169.5, -169.25, -169, -168.75, -168.5, -168.25, -168, -167.75, -167.5, -167.25, -167, -166.75, -166.5, -166.25, -166, -165.75, -165.5, -165.25, -165, -164.75, -164.5, -164.25, -164, -163.75, -163.5, -163.25, -163, -162.75, -162.5, -162.25, -162, -161.75, -161.5, -161.25, -161, -160.75, -160.5, -160.25, -160, -159.75, -159.5, -159.25, -159, -158.75, -158.5, -158.25, -158, -157.75, -157.5, -157.25, -157, -156.75, -156.5, -156.25, -156, -155.75, -155.5, -155.25, -155, -154.75, -154.5, -154.25, -154, -153.75, -153.5, -153.25, -153, -152.75, -152.5, -152.25, 88.25, 88.5, 88.75, 89, 89.25, 89.5, 89.75, 90, 90.25, 90.5, 90.75, 91, 91.25, 91.5, 91.75, 92, 92.25, 92.5, 92.75, 93, 93.25, 93.5, 93.75, 94, 94.25, 94.5, 94.75, 95, 95.25, 95.5, 95.75, 96, 96.25, 96.5, 96.75, 97, 97.25, 97.5, 97.75, 98, 98.25, 98.5, 98.75, 99, 99.25, 99.5, 99.75, 100, 100.25, 100.5, 100.75, 101, 101.25, 101.5, 101.75, 102, 102.25, 102.5, 102.75, 103, 103.25, 103.5, 103.75, 104, 104.25, 104.5, 104.75, 105, 105.25, 105.5, 105.75, 106, 106.25, 106.5, 106.75, 107, 107.25, 107.5, 107.75, 108, 108.25, 108.5, 108.75, 109, 109.25, 109.5, 109.75, 110, 110.25, 110.5, 110.75, 111, 111.25, 111.5, 111.75, 112, 112.25, 112.5, 112.75, 113, 113.25, 113.5, 113.75, 114, 114.25, 114.5, 114.75, 115, 115.25, 115.5, 115.75, 116, 116.25, 116.5, 116.75, 117, 117.25, 117.5, 117.75, 118, 118.25, 118.5, 118.75, 119, 119.25, 119.5, 119.75, 120, 120.25, 120.5, 120.75, 121, 121.25, 121.5, 121.75, 122, 122.25, 122.5, 122.75, 123, 123.25, 123.5, 123.75, 124, 124.25, 124.5, 124.75, 125, 125.25, 125.5, 125.75, 126, 126.25, 126.5, 126.75, 127, 127.25, 127.5, 127.75, 128, 128.25, 128.5, 128.75, 129, 129.25, 129.5, 129.75, 130, 130.25, 130.5, 130.75, 131, 131.25, 131.5, 131.75, 132, 132.25, 132.5, 132.75, 133, 133.25, 133.5, 133.75, 134, 134.25, 134.5, 134.75, 135, 135.25, 135.5, 135.75, 136, 136.25, 136.5, 136.75, 137, 137.25, 137.5, 137.75, 138, 138.25, 138.5, 138.75, 139, 139.25, 139.5, 139.75, 140, 140.25, 140.5, 140.75, 141, 141.25, 141.5, 141.75, 142, 142.25, 142.5, 142.75, 143, 143.25, 143.5, 143.75, 144, 144.25, 144.5, 144.75, 145, 145.25, 145.5, 145.75, 146, 146.25, 146.5, 146.75, 147, 147.25, 147.5, 147.75, 148, 148.25, 148.5, 148.75, 149, 149.25, 149.5, 149.75, 150, 150.25, 150.5, 150.75, 151, 151.25, 151.5, 151.75, 152, 152.25, 152.5, 152.75, 153, 153.25, 153.5, 153.75, 154, 154.25, 154.5, 154.75, 155, 155.25, 155.5, 155.75, 156, 156.25, 156.5, 156.75, 157, 157.25, 157.5, 157.75, 158, 158.25, 158.5, 158.75, 159, 159.25, 159.5, 159.75, 160, 160.25, 160.5, 160.75, 161, 161.25, 161.5, 161.75, 162, 162.25, 162.5, 162.75, 163, 163.25, 163.5, 163.75, 164, 164.25, 164.5, 164.75, 165, 165.25, 165.5, 165.75, 166, 166.25, 166.5, 166.75, 167, 167.25, 167.5, 167.75, 168, 168.25, 168.5, 168.75, 169, 169.25, 169.5, 169.75, 170, 170.25, 170.5, 170.75, 171, 171.25, 171.5, 171.75, 172, 172.25, 172.5, 172.75, 173, 173.25, 173.5, 173.75, 174, 174.25, 174.5, 174.75, 175, 175.25, 175.5, 175.75, 176, 176.25, 176.5, 176.75, 177, 177.25, 177.5, 177.75, 178, 178.25, 178.5, 178.75, 179, 179.25, 179.5, 179.75 ;

era5 data cf_var

Modify code on lines 96:102 to include ERA5 variables (tp, mx2t, mn2t, 2t for pr, tasmax, tasmin, tas respectively)

if variable in ['pr', 'precip', 'mtpr', 'tp']:
  cf_var = 'pr'
elif variable in ['tasmax', 'tmax', 'mx2t']:
  cf_var = 'tasmax'
elif variable in ['tasmin', 'tmin', 'mn2t']:
  cf_var = 'tasmin'
elif variable in ['tas', 't2m', '2t']:
  if time_agg == 'max':
    cf_var = 'tasmax'
  elif time_agg == 'min':
    cf_var = 'tasmin'
  else:
    cf_var = 'tas'```

icclim time bounds issues

icclim outputs an xarray Dataset with a time bounds axis. For some reason I get the following error when reading CWD index data calculated from AGCD data (i.e. /g/data/xv83/dbi599/indices/cwd_year_AGCD_v1_r005_1900-2021.nc) but not for any other metric or dataset (even other AGCD data/metrics).

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:261, in decode_cf_datetime(num_dates, units, calendar, use_cftime)
    260 try:
--> 261     dates = _decode_datetime_with_pandas(flat_num_dates, units, calendar)
    262 except (KeyError, OutOfBoundsDatetime, OverflowError):

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:217, in _decode_datetime_with_pandas(flat_num_dates, units, calendar)
    216 warnings.filterwarnings("ignore", "invalid value encountered", RuntimeWarning)
--> 217 pd.to_timedelta(flat_num_dates.min(), delta) + ref_date
    218 pd.to_timedelta(flat_num_dates.max(), delta) + ref_date

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/pandas/core/tools/timedeltas.py:148, in to_timedelta(arg, unit, errors)
    147 # ...so it must be a scalar value. Return scalar.
--> 148 return _coerce_scalar_to_timedelta_type(arg, unit=unit, errors=errors)

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/pandas/core/tools/timedeltas.py:156, in _coerce_scalar_to_timedelta_type(r, unit, errors)
    155 try:
--> 156     result = Timedelta(r, unit)
    157 except ValueError:

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/pandas/_libs/tslibs/timedeltas.pyx:1357, in pandas._libs.tslibs.timedeltas.Timedelta.__new__()

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/pandas/_libs/tslibs/timedeltas.pyx:288, in pandas._libs.tslibs.timedeltas.convert_to_timedelta64()

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/pandas/_libs/tslibs/conversion.pyx:125, in pandas._libs.tslibs.conversion.cast_from_unit()

OverflowError: Python int too large to convert to C long

During handling of the above exception, another exception occurred:

OverflowError                             Traceback (most recent call last)
File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:174, in _decode_cf_datetime_dtype(data, units, calendar, use_cftime)
    173 try:
--> 174     result = decode_cf_datetime(example_value, units, calendar, use_cftime)
    175 except Exception:

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:263, in decode_cf_datetime(num_dates, units, calendar, use_cftime)
    262 except (KeyError, OutOfBoundsDatetime, OverflowError):
--> 263     dates = _decode_datetime_with_cftime(
    264         flat_num_dates.astype(float), units, calendar
    265     )
    267     if (
    268         dates[np.nanargmin(num_dates)].year < 1678
    269         or dates[np.nanargmax(num_dates)].year >= 2262
    270     ):

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:195, in _decode_datetime_with_cftime(num_dates, units, calendar)
    193     raise ModuleNotFoundError("No module named 'cftime'")
    194 return np.asarray(
--> 195     cftime.num2date(num_dates, units, calendar, only_use_cftime_datetimes=True)
    196 )

File src/cftime/_cftime.pyx:584, in cftime._cftime.num2date()

File src/cftime/_cftime.pyx:383, in cftime._cftime.cast_to_int()

OverflowError: time values outside range of 64 bit signed integers

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Input In [24], in <cell line: 1>()
----> 1 cwd_annual_mean['AGCD'] = read_data(
      2     cwd_files['AGCD'], regrid=False, time_bounds=[start_date, end_date]
      3 )

Input In [8], in read_data(infile, regrid, time_bounds)
      1 def read_data(infile, regrid=False, time_bounds=None):
      2     """Read data and calculate annual mean.
      3     
      4     Parameters
   (...)
      9     
     10     """
---> 12     ds = xr.open_dataset(infile, decode_timedelta=False)
     13     if time_bounds:
     14         start_date, end_date = time_bounds

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/backends/api.py:495, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, backend_kwargs, *args, **kwargs)
    483 decoders = _resolve_decoders_kwargs(
    484     decode_cf,
    485     open_backend_dataset_parameters=backend.open_dataset_parameters,
   (...)
    491     decode_coords=decode_coords,
    492 )
    494 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 495 backend_ds = backend.open_dataset(
    496     filename_or_obj,
    497     drop_variables=drop_variables,
    498     **decoders,
    499     **kwargs,
    500 )
    501 ds = _dataset_from_backend_dataset(
    502     backend_ds,
    503     filename_or_obj,
   (...)
    510     **kwargs,
    511 )
    512 return ds

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:567, in NetCDF4BackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, lock, autoclose)
    565 store_entrypoint = StoreBackendEntrypoint()
    566 with close_on_error(store):
--> 567     ds = store_entrypoint.open_dataset(
    568         store,
    569         mask_and_scale=mask_and_scale,
    570         decode_times=decode_times,
    571         concat_characters=concat_characters,
    572         decode_coords=decode_coords,
    573         drop_variables=drop_variables,
    574         use_cftime=use_cftime,
    575         decode_timedelta=decode_timedelta,
    576     )
    577 return ds

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/backends/store.py:27, in StoreBackendEntrypoint.open_dataset(self, store, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta)
     24 vars, attrs = store.load()
     25 encoding = store.get_encoding()
---> 27 vars, attrs, coord_names = conventions.decode_cf_variables(
     28     vars,
     29     attrs,
     30     mask_and_scale=mask_and_scale,
     31     decode_times=decode_times,
     32     concat_characters=concat_characters,
     33     decode_coords=decode_coords,
     34     drop_variables=drop_variables,
     35     use_cftime=use_cftime,
     36     decode_timedelta=decode_timedelta,
     37 )
     39 ds = Dataset(vars, attrs=attrs)
     40 ds = ds.set_coords(coord_names.intersection(vars))

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/conventions.py:516, in decode_cf_variables(variables, attributes, concat_characters, mask_and_scale, decode_times, decode_coords, drop_variables, use_cftime, decode_timedelta)
    509     continue
    510 stack_char_dim = (
    511     concat_characters
    512     and v.dtype == "S1"
    513     and v.ndim > 0
    514     and stackable(v.dims[-1])
    515 )
--> 516 new_vars[k] = decode_cf_variable(
    517     k,
    518     v,
    519     concat_characters=concat_characters,
    520     mask_and_scale=mask_and_scale,
    521     decode_times=decode_times,
    522     stack_char_dim=stack_char_dim,
    523     use_cftime=use_cftime,
    524     decode_timedelta=decode_timedelta,
    525 )
    526 if decode_coords in [True, "coordinates", "all"]:
    527     var_attrs = new_vars[k].attrs

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/conventions.py:364, in decode_cf_variable(name, var, concat_characters, mask_and_scale, decode_times, decode_endianness, stack_char_dim, use_cftime, decode_timedelta)
    362     var = times.CFTimedeltaCoder().decode(var, name=name)
    363 if decode_times:
--> 364     var = times.CFDatetimeCoder(use_cftime=use_cftime).decode(var, name=name)
    366 dimensions, data, attributes, encoding = variables.unpack_for_decoding(var)
    367 # TODO(shoyer): convert everything below to use coders

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:673, in CFDatetimeCoder.decode(self, variable, name)
    671 units = pop_to(attrs, encoding, "units")
    672 calendar = pop_to(attrs, encoding, "calendar")
--> 673 dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
    674 transform = partial(
    675     decode_cf_datetime,
    676     units=units,
    677     calendar=calendar,
    678     use_cftime=self.use_cftime,
    679 )
    680 data = lazy_elemwise_func(data, transform, dtype)

File /g/data/xv83/dbi599/miniconda3/envs/model-eval/lib/python3.10/site-packages/xarray/coding/times.py:184, in _decode_cf_datetime_dtype(data, units, calendar, use_cftime)
    176     calendar_msg = (
    177         "the default calendar" if calendar is None else f"calendar {calendar!r}"
    178     )
    179     msg = (
    180         f"unable to decode time units {units!r} with {calendar_msg!r}. Try "
    181         "opening your dataset with decode_times=False or installing cftime "
    182         "if it is not installed."
    183     )
--> 184     raise ValueError(msg)
    185 else:
    186     dtype = getattr(result, "dtype", np.dtype("object"))

ValueError: unable to decode time units 'days since 1900-01-01 00:00:00' with "calendar 'proleptic_gregorian'". Try opening your dataset with decode_times=False or installing cftime if it is not installed.

Issue calculating CD/CW/WW/WD with CCAM ERA5

Not sure what's changed but I can't seem to calculate CD/CW/WW/WD indices in CCAM ERA5.

I think it's something to do with the lat bounds and CCAM has had some issues previously with lat/lon bounds and floating points, however I had previously calculated these indices using icclim 6.0.0 about 2 weeks ago. I checked the raw input files /g/data/xv83/mxt599/ccam_era5_evaluation_aus-10i_12km/drs_cordex/CORDEX/output/AUS-10i/CSIRO/ECMWF-ERA5/evaluation/r1i1p1f1/CSIRO-CCAM-2203/v1/day/tas and /g/data/xv83/mxt599/ccam_era5_evaluation_aus-10i_12km/drs_cordex/CORDEX/output/AUS-10i/CSIRO/ECMWF-ERA5/evaluation/r1i1p1f1/CSIRO-CCAM-2203/v1/day/pr and they haven't been changed since the start of November. Other indices calculate ok so far.

2022-12-05 09:38:09,633 Array size: (15706, 611, 928) 2022-12-05 09:38:09,633 Chunk size: Frozen({'time': (15706,), 'lat': (37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 19), 'lon': (56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 32)}) 2022-12-05 09:38:10,095 Array size: (15706, 611, 928) 2022-12-05 09:38:10,096 Chunk size: Frozen({'time': (15706,), 'lat': (37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 19), 'lon': (56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 32)}) 2022-12-05 09:38:10,096 ******************************************************************************************** 2022-12-05 09:38:10,096 * * 2022-12-05 09:38:10,096 * icclim 6.1.3 * 2022-12-05 09:38:10,096 * * 2022-12-05 09:38:10,096 * * 2022-12-05 09:38:10,096 * Sun Dec 4 22:38:10 2022 * 2022-12-05 09:38:10,096 * * 2022-12-05 09:38:10,096 * BEGIN EXECUTION * 2022-12-05 09:38:10,096 * * 2022-12-05 09:38:10,096 ******************************************************************************************** 2022-12-05 09:38:10,096 Processing: 0% Traceback (most recent call last): File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 526, in <module> main(args) File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 369, in main index = icclim.index( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/icclim/main.py", line 413, in index climate_vars_dict = read_in_files( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/icclim/models/climate_variable.py", line 223, in read_in_files return _build_in_file_dict(in_files, var_names, threshold, standard_index) File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/icclim/models/climate_variable.py", line 235, in _build_in_file_dict input_dataset = read_dataset( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/icclim/pre_processing/input_parsing.py", line 65, in read_dataset return xr.merge( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/merge.py", line 1025, in merge merge_result = merge_core( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/merge.py", line 757, in merge_core variables, out_indexes = merge_collected( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/merge.py", line 302, in merge_collected merged_vars[name] = unique_variable( File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/merge.py", line 156, in unique_variable raise MergeError( xarray.core.merge.MergeError: conflicting values for variable 'lat_bnds' on objects to be combined. You can skip this check by specifying compat='override'.

temporal aggregation

It might be useful for Climate Innovation Hub work to include a time_aggregation.py script in this repo for quickly aggregating monthly data to a custom seasons or annual timescale.

It could use xCDAT if the following issue is addressed: xCDAT/xcdat#416 (comment)

DTR/ETR/vDTR results have degK to degC temperature conversion applied

I'm not sure if this is an icclim issue or a run_icclim.py issue but I've found that DTR/vDTR/EDTR variables have 273.15 subtracted from the results. It appears that the conversion from degK to degC is being applied after the index calculation.

Output temperature units

I think the current version outputs temperature indices with units "C" whereas previously it was "°C".
For some reason this is breaking things in ILAMB. Is it possible to keep the units as "°C"? What is the standard meant to be?

  File "/home/599/bxn599/.local/lib/python3.9/site-packages/ILAMB/Variable.py", line 919, in convert
    src_unit.convert(self.data,tar_unit,inplace=True)
  File "/g/data/hh5/public/apps/miniconda3/envs/analysis3-22.10/lib/python3.9/site-packages/cf_units/__init__.py", line 1918, in convert
    raise ValueError(
ValueError: Unable to convert from 'Unit('C')' to 'Unit('°C')'.

Differences in percentiles when using --base_period with a different study period?

When calculating percentiles, icclim can use the base_period_time_range parameter to set the reference period. If base_period_time_range is not defined, the time_range is used.
I originally calculated the percentiles over the 1985-2014 time_range without base_period_time_range however if I calculate percentiles over the 1979-2021 time range with a base_period_time_range of 1985-2014, there are some very small differences. This plot shows the R95p difference for Jan 1986 between the two

Perhaps I'm misunderstanding how base_period_time_range is applied, I thought that they should be the same within 1985-2014 as the reference period is the same?

I also tested it by recalculating the percentiles over the 1985-2014 time_range with a base_period_time_range of 1985-2014 and this gave the same result as when using the 1985-2014 time_range without base_period_time_range

ERA5 data inputs

Request implementation to allow use with ERA5 data in rt52:

hourly screen temperature data in /g/data/rt52/era5/single-levels/reanalysis/2t/$year/2t_*.nc (variable name is t2m in units of K)
hourly total precipitation data in /g/data/rt52/era5/single-levels/reanalysis/mtpr/$year/mtpr_*.nc (variable name is mtpr in kg m**-2 s**-1units of )

I assume icclim does not handle hourly dat, so an additional argument can be passed to run_icclim.py, which tells xr.dataarray how to preprocess the hourly data first such as resample('1D').mean() etc.

Cannot run on GCMs NorESM2-MM and CMCC-ESM2

I've started running icclim on the GCMs being downscaled and have run into issues with NorESM2-MM and CMCC-ESM2. Looks like its the same error for the two models. I have successfully run icclim on the other GCMs and will be adding these to ia39 test-dir shortly.

CMCC-ESM2:

Traceback (most recent call last):
  File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 536, in <module>
    main(args)
  File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 365, in main
    ds, cf_var = read_data(
  File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 283, in read_data
    ds = xr.open_mfdataset(infiles, chunks='auto', mask_and_scale=True)
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 996, in open_mfdataset
    datasets = [open_(p, **open_kwargs) for p in paths]
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 996, in <listcomp>
    datasets = [open_(p, **open_kwargs) for p in paths]
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 545, in open_dataset
    ds = _dataset_from_backend_dataset(
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 357, in _dataset_from_backend_dataset
    ds = _chunk_ds(
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 325, in _chunk_ds
    var_chunks = _get_chunk(var, chunks)
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/dataset.py", line 220, in _get_chunk
    chunk_shape = da.core.normalize_chunks(
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py", line 3073, in normalize_chunks
    chunks = auto_chunks(chunks, shape, limit, dtype, previous_chunks)
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py", line 3169, in auto_chunks
    raise NotImplementedError(
NotImplementedError: Can not use auto rechunking with object dtype. We are unable to estimate the size in bytes of object data

NorESM2-MM:

Traceback (most recent call last):
  File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 536, in <module>
    main(args)
  File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 365, in main
    ds, cf_var = read_data(
  File "/g/data/xv83/bxn599/ACS/icclim/run_icclim.py", line 283, in read_data
    ds = xr.open_mfdataset(infiles, chunks='auto', mask_and_scale=True)
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 996, in open_mfdataset
    datasets = [open_(p, **open_kwargs) for p in paths]
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 996, in <listcomp>
    datasets = [open_(p, **open_kwargs) for p in paths]
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 545, in open_dataset
    ds = _dataset_from_backend_dataset(
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 357, in _dataset_from_backend_dataset
    ds = _chunk_ds(
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/backends/api.py", line 325, in _chunk_ds
    var_chunks = _get_chunk(var, chunks)
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/xarray/core/dataset.py", line 220, in _get_chunk
    chunk_shape = da.core.normalize_chunks(
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py", line 3073, in normalize_chunks
    chunks = auto_chunks(chunks, shape, limit, dtype, previous_chunks)
  File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/dask/array/core.py", line 3169, in auto_chunks
    raise NotImplementedError(
NotImplementedError: Can not use auto rechunking with object dtype. We are unable to estimate the size in bytes of object data

TNx as TNn

We (Emma) have spotted a curious case where running,

# Example: TNx: Maximum daily minimum temperature
input_files=/g/data/xv83/agcd-csiro/tmin/daily/tmin_AGCD-CSIRO_r005_20200101-20201231_daily.nc
/g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python /g/data/tp28/dev/chs548/indices/run_icclim.py --slice_mode month --verbose --input_files ${input_files} --variable tmin TNx icclim.out.nc

Produces a netcdf file that contains,

  variables:
    float TNx(time,lat,lon) ;
      TNx:_FillValue = NaNf ;
      TNx:standard_name = "minimum_air_temperature" ;
      TNx:long_name = "Minimum of daily minimum air temperature for each month." ;
      string TNx:units = "°C" ;

It should be MAXIMUM of daily minimum air temperature. We also found that the resulting values "TNx" are equivalent to "TNn".

My test case here: gadi: /scratch/public/chs548/test/

Output for R75pTOT/R95pTOT/R99pTOT is empty for CCAM/BARPA/CMIP6 models

icclim isn't outputting data for R75pTOT/R95pTOT/R99pTOT for CCAM/BARPA or CMIP6 models. It worked fine for AGCDv1. I am wondering if this is caused by unit differences? AGCDv1 is mm and CCAM/CMIP6/BARPA is kg m-2 s-1.

This plot is from /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/AUS-15/BOM/ECMWF-ERA5/evaluation/none/BOM-BARPA-R/v1/climdex/R75pTOT/R75pTOT_AUS-15_ECMWF-ERA5_evaluation_v1_month_19850101-20141231.nc

This plot is from /g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/indices/AUS-r005/none/BOM-AGCD/historical/none/none/none/climdex/R99pTOT/R99pTOT_AUS-r005_BOM-AGCD_historical_none_month_19850101-20141231.nc

Output percentile indices with (time, lat, lon) order

Is it possible for percentile indices to be output with time as the left-most dimension? Currently output is (lat, lon, time) which doesn't seem to play nicely with ILAMB. It only occurs for percentile indices, other indices are output with time as the left-most dimension.

icclim lat/lon subset

@DamienIrving .

When running

/g/data/xv83/dbi599/miniconda3/envs/icclim/bin/python /g/data/tp28/dev/chs548/indices/run_icclim.py --slice_mode month --verbose --start_date 1980-01-01 --end_date 2020-12-31 --input_files /g/data/xv83/bxn599/test-data/CORDEX-CMIP6/output//AUS-10i/CSIRO/ECMWF-ERA5/evaluation/r1i1p1f1/CSIRO-CCAM-2203/v1/day/tasmax/tasmax_AUS-10i_ECMWF-ERA5_evaluation_r1i1p1f1_CSIRO-CCAM-2203_v1_day_19800101-19801231.nc --variable tasmax tx out.nc

led to,

Traceback (most recent call last): File "/g/data/tp28/dev/chs548/indices/run_icclim.py", line 391, in <module> main(args) File "/g/data/tp28/dev/chs548/indices/run_icclim.py", line 229, in main ds, cf_var = read_data( File "/g/data/tp28/dev/chs548/indices/run_icclim.py", line 167, in read_data ds = subset_bbox(ds, start_date=start_date, end_date=end_date, lat_bnds=lat_bnds, lon_bnds=lon_bnds) File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/clisops/core/subset.py", line 270, in func_checker return func(*args, **kwargs) File "/g/data/xv83/dbi599/miniconda3/envs/icclim/lib/python3.10/site-packages/clisops/core/subset.py", line 1211, in subset_bbox da = da.sel({lon: slice(*lon_bnds)}) TypeError: iteration over a 0-d array

Same when running on BARPA data.

Not obvious to me - possible bug in subset.py, since the icclim.py is correctly passing lat_bnds=None and lon_bnds=None to subset_bbox function.

ausclimateservice / indices Goto Github PK

indices's Introduction

indices

Python environment

Usage: run_icclim.py

Example 1: Simple indices

Example 2: Sub-daily input data

Example 3: Multi-variate indices

Example 4: "Along the time axis" indices

Usage: run_xclim.py

Example 1: Simple indices

Example 2: Indices involving daily mean temperature

Common configurations

indices's People

Contributors

Stargazers

Watchers

Forkers

indices's Issues

Recommend Projects

Recommend Topics

Recommend Org

Usage: `run_icclim.py`

Usage: `run_xclim.py`