Rechunk DataArray to calculate 90% quantile over over chunked time dimension

Viewed 41

I try to calculate the 90 percentile of a variable over a period of 16 years. The data is stored in netCDF files (where 1 month is stored in 1 file --> 12files/year * 16years).

I pre-processed the data and took the daily_max and monthly mean of the variable of interested. So bottom line the folder consists of 192 files that contain each one value (the monthly mean of the daily max).

The data was opened using following command:

ds = xr.open_mfdataset(f"{folderdir}/*.nc", chunks={"time":1})

Trying to calculate the quantile (from some data variable, which was extracted from the ds: data_variable = ds["data_variable"]) with following code:

q90 = data_varaible.qunatile(0.95, "time"), yields follwing error message:

ValueError: dimension time on 0th function argument to apply_ufunc with dask='parallelized' consists of multiple chunks, but is also a core dimension. To fix, either rechunk into a single dask array chunk along this dimension, i.e., .chunk(dict(time=-1)), or pass allow_rechunk=True in dask_gufunc_kwargs but beware that this may significantly increase memory usage.

I tried to rechunk, as explained in the error message by apply: data_variable.chunk(dict(time=-1).quantile(0.95,'time'), with no success (got the exact same error. Further I tired to rechunk in the following way: data_variable.chunk({'time':1})), which was also not successful.

Printing out the data.variable.chunk(), actually shows that the chunk size in time dimension is supposed to be 1, so i don't understand where I made a mistake.

ps: I didn't try allow_rechunk=True in dask_gufunc_kwargs, since I don't know where to pass that argument.

Thanks for the help,

Max

ps: Printing out the data_variable yields, (to be clear, some_variable (see above) is 'wsgsmax' here):

<xarray.DataArray 'wsgsmax' (time: 132, y: 853, x: 789)>
dask.array<concatenate, shape=(132, 853, 789), dtype=float32, chunksize=(1, 853, 789), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 1995-01-16T12:00:00 ... 2005-12-16T12:00:00
    lon      (y, x) float32 dask.array<chunksize=(853, 789), meta=np.ndarray>
    lat      (y, x) float32 dask.array<chunksize=(853, 789), meta=np.ndarray>
  * x        (x) float32 0.0 2.5e+03 5e+03 ... 1.965e+06 1.968e+06 1.97e+06
  * y        (y) float32 0.0 2.5e+03 5e+03 ... 2.125e+06 2.128e+06 2.13e+06
    height   float32 10.0
Attributes:
    standard_name:  wind_speed_of_gust
    long_name:      Maximum Near Surface Wind Speed Of Gust
    units:          m s-1
    grid_mapping:   Lambert_Conformal
    cell_methods:   time: maximum
    FA_name:        CLSRAFALES.POS
    par:            228
    lvt:            105
    lev:            10
    tri:            2
1 Answers

chunk({"time": 1} will produce as many chunks as there are time steps. Each chunk will have a size of 1.

Printing out the data.variable.chunk(), actually shows that the chunk size in time dimension is supposed to be 1, so i don't understand where I made a mistake.

To compute percentiles dask needs to load the full timeseries into memory so it forbids chunking over "time" dimension. So what you want is either chunk({"time": len(ds.time)} or to use directly the shorthand chunk({"time": -1}.

I don't understand why data_variable.chunk(dict(time=-1).quantile(0.95,'time') would not work though.

Related