Code Sample

```python def cut_dataset(ds2cut, varList = ['Temp' 'S' 'Eta' 'U' 'V' 'W'], lonRange = [-180, 180], latRange = [-90, 90], depthRange = [0, float("inf")], timeRange = ['2007-09-01T00', '2008-08-31T18'], timeFreq = '1D', sampMethod = 'mean', interpC = True, saveNetCDF = False): """ Cut the dataset """

# Copy dataset
ds = ds2cut.copy(deep=True)

# Choose variables
varList_tmp = list(varList)
for varName in ds.variables:
    if all(x != 'time' for x in ds[varName].dims) or (varName=='time'):
            varList_tmp.append(varName)
toDrop = list(set(ds.variables)-set(varList_tmp))
ds = ds.drop(toDrop)

# Cut dataset
ds = ds.sel(time = slice(min(timeRange),  max(timeRange)),
            Xp1  = slice(min(lonRange),   max(lonRange)),
            Yp1  = slice(min(latRange),   max(latRange)),
            Zp1  = slice(min(depthRange), max(depthRange)))
ds = ds.sel(X    = slice(min(ds['Xp1'].values), max(ds['Xp1'].values)),
            Y    = slice(min(ds['Yp1'].values), max(ds['Yp1'].values)),
            Z    = slice(min(ds['Zp1'].values), max(ds['Zp1'].values)),
            Zu   = slice(min(ds['Zp1'].values), max(ds['Zp1'].values)),
            Zl   = slice(min(ds['Zp1'].values), max(ds['Zp1'].values)))

# Resample
if sampMethod=='snapshot':
    ds = ds.resample(time=timeFreq).first(skipna=False)
elif sampMethod=='mean':
    ds = ds.resample(time=timeFreq).mean()

# Create grid
grid = xgcm.Grid(ds,periodic=False)

# Interpolate
if interpC:
    for varName in varList:
        for dim in ds[varName].dims:
            if len(dim)>1 and dim!='time':
                ds[varName] = grid.interp(ds[varName], axis=dim[0])

# Remove useless variables
allDims = []
for varName in varList:
    for dim in ds[varName].dims:
        allDims.append(dim)
toDrop = []
for varName in ds.variables:
    if len(list(set(ds[varName].dims)-set(allDims)))>0:
        toDrop.append(varName)
ds = ds.drop(toDrop)

# Save to NetCDF
if saveNetCDF: ds.to_netcdf(saveNetCDF)

return ds, grid

3D test

ds_cut, grid_cut = cut_dataset(ds, varList = ['Eta'], latRange = [65, 65.5], depthRange = [0, 2], timeRange = ['2007-11-15T00', '2007-11-16T00'], timeFreq = '1D', sampMethod = 'mean', interpC = False, saveNetCDF = '3Dvariable.nc')

4D test

ds_cut, grid_cut = cut_dataset(ds, varList = ['Temp'], lonRange = [-30, -29.5], latRange = [65, 65.5], depthRange = [0, 2], timeRange = ['2007-11-15T00', '2007-11-16T00'], timeFreq = '1D', sampMethod = 'mean', interpC = False, saveNetCDF = '4Dvariable.nc') ```

Problem description

I'm working with a big dataset. However, most of the time I only need a small subset of data. My idea was to open and concatenate everything with open_mfdataset, and then extract subsets of data using the indexing routines. This approach works very good when I extract 3D variables (just lon, lat, and time), but it fails when I try to extract 4D variables (lon, lat, time, and depth). It doesn't actually fail, but to_netcdf takes forever. When I open a smaller dataset since the very beginning (let's say just November), then I'm also able to extract 4D variables. When I load the sub-dataset after using the indexing routines, does xarray need to read the whole original 4D variable? If yes, then I should probably change my approach and I should open subset of data since the very beginning. If no, am I doing something wrong?

Output of xr.show_versions()