Create one NetCDF file from multiple NetCDF files using Python [closed]

Question

I am downloading climate data in NetCDF format. For each variable (e.g. 'precipitation'), I need to merge 9 NetCDFs, each belonging to an unique climate model. Each NetCDF has the same size (time, lat, lon).

How can I merge 9 3D NetCDFs into one 4D NetCDF?

Ultimately, I want to calculate cumulative precipitation per month. Here's my current code:

variables = ['pr']         
scenarios = ['historical', 'ssp245']        #options ['historical', 'ssp126', 'ssp245', 'ssp370', 'ssp585']
models = ['UKESM1-0-LL', 'MRI-ESM2-0', 'MIROC6', 'MIROC-ES2L', 'IPSL-CM6A-LR',
         'GFDL-ESM4', 'FGOALS-g3', 'CNRM-ESM2-1', 'CanESM5']


save_folder = processing_fn / 'local_climate_assessment' / f'{variable}' / 'output'
if not os.path.exists(save_folder):
    os.makedirs(save_folder)

netcdfs = []

# Create one netcdf per model by merging annual netcdfs
for variable in variables:
    for scenario in scenarios:
        for model in models:
    

            source = processing_fn / 'local_climate_assessment' / f'{variable}' / f'{scenario}' / f'{model}'
            netcdf_fn = save_folder / f'{variable}_{scenario}_{model}.nc'
            
            if not os.path.exists(netcdf_fn):
        
                gdf_model = xr.open_mfdataset(str(source / '*.nc'), combine = 'nested', concat_dim="time", use_cftime=True)
                # rename_dict = {variable, f'{variable}_{scenario}_{model}'}
                # gdf_model.rename(rename_dict, inplace = True)
                gdf_model.to_netcdf(netcdf_fn)
                print(gdf_model.attrs['cmip6_source_id'])
                netcdfs.append(gdf_model)
                
            else:
                gdf_model = xr.open_mfdataset(netcdf_fn)
                netcdfs.append(gdf_model)

# Create one netcdf per variable by merging models
ds = xr.combine_nested(netcdfs, concat_dim = "time")
print(ds)
Out[33]: 
<xarray.Dataset>
Dimensions:  (time: 246095, lat: 47, lon: 50)
Coordinates:
  * time     (time) object 1981-01-01 12:00:00 ... 2060-12-31 12:00:00
  * lat      (lat) float64 31.62 31.88 32.12 32.38 ... 42.38 42.62 42.88 43.12
  * lon      (lon) float64 234.6 234.9 235.1 235.4 ... 246.1 246.4 246.6 246.9
Data variables:
    pr       (time, lat, lon) float32 dask.array<chunksize=(360, 47, 50), meta=np.ndarray>

The above code works, but I'm creating one big 3D NetCDF instead of a 4D still containing the climate model names. The code below results in the following error:

a = ds.resample(time = 'M').sum()
ValueError: index must be monotonic for resampling

How to create a 4D NetCDF with model names included, and resample to create monthly sum values?

Answer 1

The trick is to first add the new (4th) dimension to each individual dataset, then you can easily merge them. Here I do this for just 2 datasets, but I think you can easily expand it to 9.

import xarray as xr

# Create two example 3D-datasets.
ds1 = xr.Dataset(None, {"time": None, "x": None, "y": None})
ds2 = xr.Dataset(None, {"time": None, "x": None, "y": None})

# Make a list of the datasets and make a list with their respective names.
model_dss = [ds1, ds2]
model_names = ["model1", "model2"]

# Add a new dimension to each dataset.
expanded_dss = [ds.expand_dims({"model": [model_name]}) for ds, model_name in zip(model_dss, model_names)]

# Merge the two 4D datasets together.
ds_4D = xr.merge(expanded_dss)