Functions for retrieving the time axis from multiple netCDF data files that do not have the same time units

mfdataset_time.py

"""Functions for retrieving the time axis from multiple
netCDF data files that do not have the same time units"""
import netCDF4
import glob
import numpy as np


def read_mfdataset_time(pattern,timevar='time',calendar='proleptic_gregorian'):
    """Get time axis from multi-file netCDF data using np.append"""
    files = sorted(glob.glob(pattern))
    for i,fname in enumerate(files):
        with netCDF4.Dataset(fname) as ds:
            if i == 0:
                singleshape = ds.variables[timevar].shape
                time = np.zeros([0]*len(singleshape),object)
            time = np.append(time,netCDF4.num2date(
                ds.variables[timevar][:],
                ds.variables[timevar].units,calendar),axis=0)
    return time


def read_mfdataset_time_same_length(pattern,timevar='time',calendar='proleptic_gregorian'):
    """Get time axis from multi-file netCDF data assuming same length for all files"""
    files = sorted(glob.glob(pattern))
    for i,fname in enumerate(files):
        with netCDF4.Dataset(fname) as ds:
            if i == 0:
                singleshape = ds.variables[timevar].shape
                nt = singleshape[0]
                targetshape = [len(files)*nt]
                if len(singleshape) > 1:
                    targetshape.append(singleshape[1:])
                time = np.zeros(targetshape,object)
            ii = i*nt
            time[ii:ii+nt] = netCDF4.num2date(
                ds.variables[timevar][:],
                ds.variables[timevar].units,calendar)
    return time


def read_mfdataset_time_get_shapes_first(pattern,timevar='time',calendar='proleptic_gregorian'):
    """Get time axis from multi-file netCDF data reading the shapes from all files first"""
    files = sorted(glob.glob(pattern))

    nt = np.zeros(len(files))
    for i,fname in enumerate(files):
        with netCDF4.Dataset(fname) as ds:
            shape = np.array(ds.variables[timevar].shape)
            nt[i] = shape[0]
    targetshape = [np.sum(nt)]
    if len(shape) > 1: 
        targetshape.append(shape[1:])

    time = np.zeros(targetshape,object)
    ii = 0
    for i,fname in enumerate(files):
        with netCDF4.Dataset(fname) as ds:
            time[ii:ii+nt[i]] = netCDF4.num2date(
                ds.variables[timevar][:],
                ds.variables[timevar].units,calendar)
            ii += nt[i]
    return time


def demonstrate_problem(pattern):
    """Demonstrate the problem with multi-file datasets that 
    have different reference dates for their time axis.
    
    Parameters
    ----------
    pattern : str
        File search pattern to netCDF files
    """
    for fname in sorted(glob.glob(pattern)):
        with netCDF4.Dataset(fname) as ds:
            print(ds.variables['time'].units)

    with netCDF4.MFDataset(pattern) as ds:
        print(ds.variables['time'].units)



if __name__ == '__main__':

    pattern = '*.nc'

    time = read_mfdataset_time(pattern,'time')