def mean_by_batch(da, batch_size, axis=0):
    """
    Loop over batches indexed in axis in an xarray DataArray to take the grand mean of the array in a memory-
    efficient way.

    :param da: xarray DataArray
    :param batch_size: int: number of samples to load and mean at a time
    :param axis: int: axis along which to index batches
    :return: float: the mean of the array
    """
    size = da.shape[axis]
    batches = list(range(0, size, batch_size))
    dim = da.dims[axis]
    total = 0.0
    for b in batches:
        total += da.isel(**{dim: slice(b, min(b+batch_size, size))}).values.sum()
    return total / da.size