cerfacs-globc / icclim

icclim: Python library for climate indices and climate indicators calculation.
https://icclim.readthedocs.io/en/latest/
Apache License 2.0
80 stars 32 forks source link

BUG: time coordinate not properly set to the middle of the sampling period #273

Closed bzah closed 9 months ago

bzah commented 1 year ago

Description

Initially reported by email.

I am using your library to calculate the SPI6 climate index but I don't understand the results it gives me, I don't know if there is a bug. I send the netcdf that I'm using and this is the code that I'm using from you:

import icclim
dataset = icclim.index(
    index_name="spi6",
    in_files=dataset,
    var_name="pr",
    base_period_time_range=baseline_period_time_range,
).load()

and the result looks like this:

<xarray.Dataset>
Dimensions:      (lat: 192, lon: 288, time: 12, bounds: 2)
Coordinates:
  * lat          (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0
  * lon          (lon) float64 -178.8 -177.5 -176.2 -175.0 ... 177.5 178.8 180.0
    height       float64 2.0
  * time         (time) datetime64[ns] 2015-07-02 ... 2015-12-16
  * bounds       (bounds) int64 0 1
Data variables:
    SPI6         (time, lat, lon) float64 nan nan nan nan ... nan nan nan nan
    time_bounds  (time, bounds) datetime64[ns] 2015-01-01 ... 2015-12-31
Attributes:
    title:        spi
    references:   ATBD of the ECA&D indices calculation (https://knmi-ecad-as./..
    institution:  Climate impact portal (https://climate4impact.eu/)
    history:      [2023-05-25 14:06:10] spi: SPI(pr=pr, pr_cal=pr, freq='MS',...
    source:      
    Conventions:  CF-1.6 

The time dimension is supposed to start in january.

The bug is due to the post_processing step that is meant to add time_bounds but also resample the result if necessary. This step should probably be avoided when computing SPI6 or SPI3. Investigation needed to see if other indices are impacted.

bzah commented 1 year ago

Workaround until the issue is fixed, using xclim:

import xclim
import xarray as xr

ds = xr.open_dataset("pr.nc")
base_pr = ds.pr.sel(time=slice("2015-01-01","2015-06-31"))
xclim.atmos.standardized_precipitation_index(
                pr=ds.pr,
                pr_cal=base_pr,
                freq="MS", 
                window=6, 
                dist="gamma", 
                method="APP"
)
bzah commented 1 year ago

I'm working on a fix. @pagecp do you think this issue also apply to other indices ?

pagecp commented 1 year ago

I'm working on a fix. @pagecp do you think this issue also apply to other indices ?

I need to go through all the indices, to be sure.