Saving derived fields - x4c: Xarray for CESM

Writing a derived field back to disk is easy to get wrong by hand.

import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import nc_time_axis  # registers the cftime axis converter for matplotlib
import x4c

x4c.set_style('journal')

# after set_style: it resets rcParams from matplotlibrc defaults, which includes the
# backend, so assert the inline backend last or figures are never captured
%matplotlib inline

# The tutorial runs against a reduced copy of a real CESM case. It is published as a
# GitHub Release asset rather than committed, so the first call downloads it into
# ~/.cache/x4c (override with $X4C_CACHE_DIR) and later calls reuse it. Set
# $X4C_SAMPLE_DIR to point at a copy you already have.
case_dir = x4c.fetch_sample_data(case='cesm1', verbose=False)
casename = os.path.basename(case_dir)
print('x4c', x4c.__version__)

x4c 2026.6.11

case = x4c.Timeseries(case_dir, grid_dict={'atm': 'ne16np4', 'ocn': 'g16'},
                      cesm_ver=1)
case.load('TS', timespan=(1, 10), verbose=False)

# EOFs need a regular grid, so regrid off the spectral-element grid first
ts = case.ds['TS'].x.regrid(2, 2).x.da
print(dict(ts.sizes))

>>> case.root_dir: /glade/u/home/fengzhu/.cache/x4c/sample_data/cesm1_sample_data/b.e13.B1850C5.ne16_g16.icesm131_d18O_fixer.Miocene.3xCO2.005
>>> case.path_pattern: comp/proc/tseries/*/casename.hstr.vn.timespan.nc
>>> case.grid_dict: {'atm': 'ne16np4', 'ocn': 'g16', 'lnd': 'ne16np4', 'rof': 'ne16np4', 'ice': 'g16'}
>>> case.casename: b.e13.B1850C5.ne16_g16.icesm131_d18O_fixer.Miocene.3xCO2.005
>>> case.paths["atm"]["cam.h0"] created
>>> case.paths["ocn"]["pop.h"] created
>>> case.paths["lnd"]["clm2.h0"] created
>>> case.paths["ice"]["cice.h"] created
>>> case.vns["atm"]["cam.h0"] created
>>> case.vns["ocn"]["pop.h"] created
>>> case.vns["lnd"]["clm2.h0"] created
>>> case.vns["ice"]["cice.h"] created

{'time': 120, 'lat': 90, 'lon': 180}

Saving derived fields¶

Use .x.to_netcdf(), not the bare xarray method. x4c carries the grid metadata (gw, lat, lon, dz) as .attrs, and those are DataArrays — netCDF attributes have to be scalars or strings, so a plain to_netcdf fails.

.x.to_netcdf() strips them from a copy, so the object you still hold keeps its weights and stays usable afterwards.

gmst = case.calc('GMST ~ TS:ann:gm', timespan=(1, 10), verbose=False)
out_path = 'gmst_ann.nc'

gmst.x.to_netcdf(out_path)
print('wrote', out_path, f'({os.path.getsize(out_path) / 1024:.1f} KB)')

# the in-memory object is untouched
print('still has gw? ', 'gw' in case.ds['TS'].x.da.attrs)

>>> Timespan: [0001-01-01 00:00:00, 0010-12-01 00:00:00]

wrote gmst_ann.nc (9.9 KB)
still has gw?  True

back = xr.open_dataset(out_path)
print(back)
back.close()
os.remove(out_path)

<xarray.Dataset> Size: 160B
Dimensions:  (time: 10)
Coordinates:
  * time     (time) object 80B 0001-12-31 00:00:00 ... 0010-12-31 00:00:00
Data variables:
    GMST     (time) float64 80B ...

A plain to_netcdf on a field that still carries the grid attrs raises — worth seeing once so the error is recognisable.

da = case.ds['TS'].x.da
try:
    da.to_netcdf('should_fail.nc')
except TypeError as e:
    print('TypeError:', str(e)[:120], '...')
finally:
    if os.path.exists('should_fail.nc'):
        os.remove('should_fail.nc')

TypeError: Invalid value for attr 'gw': <xarray.DataArray 'area' (ncol: 13826)> Size: 111kB
dask.array<where, shape=(13826,), dtype ...