using chunks='auto' inside open_mfdataset call with a cftime coordinate

[tlsw231]

I am trying to process some massively large climate timeseries (~7 TB) and would like to use chunks='auto' inside the open_mfdataset call. In the past, including this one argument in open_mfdataset has proven to be absolutely essential while trying to process large number of files. Without that argument, dask wants to load too much of the dataset into memory which leads to failure.

My problem is that the present dataset has a cftime time coordinate, which according to this thread prevents automatic chunking. I get the same error as in that thread:

NotImplementedError: Can not use auto rechunking with object dtype. We are unable to estimate the size in bytes of object data

Is there any way out?

[spencerkclark]

I think this suggests you have at least one data variable that contains cftime.datetime objects in your Dataset (likely in addition to a time coordinate, which will automatically be loaded into memory and not have this issue). A couple options would be:

1. Drop these data variables within your open_mfdataset call (if you don't need them and you know their names):

ds = xr.open_mfdataset(..., chunks="auto", drop_variables=["time_data_var1", "time_data_var2"])

2. Or hold off on decoding the times until after you open the dataset:

ds = xr.open_mfdataset(..., chunks="auto", decode_cf=False)
ds = xr.decode_cf(ds)

This gets a little more complicated if the times in your multi-file dataset are encoded with different units. In that case you need to provide a preprocess function that ensures all the times are encoded with the same units. The code below should do that in a general way:

import cftime

def rebase_times_numpy(arr, units, calendar, target_units):
    times = cftime.num2date(arr, units, calendar)
    return cftime.date2num(times, target_units, calendar)

def rebase_times(da, units, calendar, target_units="microseconds since 1970-01-01"):
    kwargs = {"units": units, "calendar": calendar, "target_units": target_units}
    result = xr.apply_ufunc(rebase_times_numpy, da, dask="parallelized", kwargs=kwargs)
    return result.assign_attrs(units=target_units, calendar=calendar)

def _harmonize_time_like_dataarrays(d):
    time_like_dataarrays = {}
    for name, da in d.items():
        if "units" in da.attrs and "since" in da.attrs["units"]:
            units = da.attrs["units"]
            calendar = da.attrs.get("calendar", "standard")
            time_like_dataarrays[name] = rebase_times(da, units, calendar)
    return time_like_dataarrays

def preprocess(ds):
    time_like_data_vars = _harmonize_time_like_dataarrays(ds.data_vars)
    time_like_coords = _harmonize_time_like_dataarrays(ds.coords)
    return ds.assign(time_like_data_vars).assign_coords(time_like_coords)

ds = xr.open_mfdataset(..., chunks="auto", decode_cf=False, preprocess=preprocess)
ds = xr.decode_cf(ds)

I'm just guessing here, though -- it might be helpful to see the output of ds.info() for one of the files you are loading.

[tlsw231]

Thanks!

Option 1 was not possible as the timeseries has basically one time-dependent variable that I need.
Option 2 might be what I was looking for. Looks promising even if my job has a ways to go (might need to request more memory and rerun).

.