Python Climate Model Visualization — Deep Dive

Technical foundation

Climate model visualization requires handling multi-dimensional datasets (time × latitude × longitude × pressure level × model × scenario), applying correct map projections, using scientifically appropriate colormaps, and creating graphics that accurately represent uncertainty. This deep dive covers the code patterns, design decisions, and technical details for creating research-grade and communication-quality climate figures.

Accessing CMIP6 model output

The Coupled Model Intercomparison Project Phase 6 (CMIP6) is the current generation of coordinated climate model experiments. Data is distributed via the Earth System Grid Federation (ESGF):

import xarray as xr
import intake

# Using intake-esm catalog for CMIP6 data discovery
catalog = intake.open_esm_datastore(
    "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
)

# Search for surface temperature projections
results = catalog.search(
    experiment_id=["ssp245", "ssp585"],
    variable_id="tas",  # Near-surface air temperature
    table_id="Amon",    # Monthly atmospheric data
    member_id="r1i1p1f1",
    source_id=["CESM2", "UKESM1-0-LL", "MPI-ESM1-2-HR", "GFDL-ESM4"],
)

# Load into xarray datasets (lazy via Dask)
datasets = results.to_dataset_dict(zarr_kwargs={"consolidated": True})

Regridding to common resolution

Different models use different grids. Before computing multi-model statistics, regrid to a common resolution:

import xesmf as xe

# Target grid: 1° × 1° regular
target_grid = xr.Dataset({
    "lat": (["lat"], np.arange(-89.5, 90, 1.0)),
    "lon": (["lon"], np.arange(0.5, 360, 1.0)),
})

def regrid_dataset(ds: xr.Dataset, target: xr.Dataset) -> xr.Dataset:
    """Regrid a climate dataset to a common resolution."""
    regridder = xe.Regridder(ds, target, "bilinear", periodic=True)
    return regridder(ds)

# Regrid all models
regridded = {name: regrid_dataset(ds, target_grid) for name, ds in datasets.items()}

Spatial maps with Cartopy

Temperature anomaly map

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cmocean

def plot_warming_map(
    anomaly: xr.DataArray,
    title: str,
    projection=ccrs.Robinson(),
    vmin: float = -2,
    vmax: float = 8,
    cmap: str = "cmocean.thermal",
):
    """Publication-quality warming anomaly map."""
    fig, ax = plt.subplots(
        figsize=(14, 8),
        subplot_kw={"projection": projection},
    )

    # Use diverging normalization centered at 0
    norm = mcolors.TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)

    im = anomaly.plot.pcolormesh(
        ax=ax,
        transform=ccrs.PlateCarree(),
        cmap=cmocean.cm.balance,
        norm=norm,
        add_colorbar=False,
    )

    # Colorbar with proper formatting
    cbar = plt.colorbar(
        im, ax=ax, orientation="horizontal",
        fraction=0.046, pad=0.06, extend="both",
    )
    cbar.set_label("Temperature change (°C)", fontsize=12)
    cbar.ax.tick_params(labelsize=10)

    # Map features
    ax.coastlines(linewidth=0.6, color="gray")
    ax.add_feature(cfeature.BORDERS, linewidth=0.3, color="gray")
    ax.set_global()
    ax.set_title(title, fontsize=14, fontweight="bold", pad=15)

    plt.tight_layout()
    return fig

# Compute anomaly: 2070-2099 minus 1981-2010
future = ensemble_mean.sel(time=slice("2070", "2099")).mean("time")
baseline = ensemble_mean.sel(time=slice("1981", "2010")).mean("time")
warming = future - baseline

fig = plot_warming_map(warming, "Projected Warming: 2070–2099 vs 1981–2010 (SSP2-4.5)")
fig.savefig("warming_map_ssp245.png", dpi=300, bbox_inches="tight")

Stippling for statistical significance

Climate publications add stippling (dots) to show where changes are statistically significant:

from scipy import stats

def significance_mask(
    model_changes: xr.DataArray,  # (model, lat, lon)
    threshold: float = 0.05,
) -> xr.DataArray:
    """Test if multi-model mean change is significantly different from zero."""
    t_stat, p_value = stats.ttest_1samp(model_changes.values, 0, axis=0)
    significant = p_value < threshold
    return xr.DataArray(significant, coords=model_changes.coords[1:])

# Add stippling overlay
sig = significance_mask(all_model_changes)
ax.contourf(
    sig.lon, sig.lat, sig.values,
    levels=[0.5, 1.5],
    hatches=[".."],
    colors="none",
    transform=ccrs.PlateCarree(),
)

Time series with ensemble spread

def plot_ensemble_timeseries(
    scenarios: dict[str, xr.DataArray],  # scenario_name -> (model, time)
    reference_period: tuple[str, str] = ("1981", "2010"),
    title: str = "Global Mean Temperature Projections",
):
    """Multi-scenario time series with ensemble spread."""
    fig, ax = plt.subplots(figsize=(12, 6))

    scenario_colors = {
        "ssp126": "#1b9e77",
        "ssp245": "#d95f02",
        "ssp370": "#e7298a",
        "ssp585": "#7570b3",
    }

    for scenario, data in scenarios.items():
        # Compute anomaly relative to reference period
        baseline = data.sel(time=slice(*reference_period)).mean("time")
        anomaly = data - baseline

        # Area-weighted global mean
        weights = np.cos(np.deg2rad(anomaly.lat))
        global_mean = anomaly.weighted(weights).mean(["lat", "lon"])

        # Multi-model statistics
        median = global_mean.median("model")
        p10 = global_mean.quantile(0.1, "model")
        p90 = global_mean.quantile(0.9, "model")

        color = scenario_colors.get(scenario, "gray")
        time = median.time.dt.year

        ax.plot(time, median, color=color, linewidth=2, label=scenario.upper())
        ax.fill_between(time, p10, p90, color=color, alpha=0.2)

    ax.axhline(y=1.5, color="red", linestyle="--", alpha=0.5, label="1.5°C threshold")
    ax.axhline(y=2.0, color="darkred", linestyle="--", alpha=0.5, label="2.0°C threshold")

    ax.set_xlabel("Year", fontsize=12)
    ax.set_ylabel("Temperature anomaly (°C)", fontsize=12)
    ax.set_title(title, fontsize=14, fontweight="bold")
    ax.legend(loc="upper left", fontsize=10)
    ax.grid(alpha=0.3)

    plt.tight_layout()
    return fig

Warming stripes

def warming_stripes(
    annual_temp: xr.DataArray,
    reference_period: tuple[str, str] = ("1981", "2010"),
    cmap: str = "RdBu_r",
):
    """Generate Hawkins-style warming stripes."""
    baseline = annual_temp.sel(time=slice(*reference_period)).mean()
    anomaly = annual_temp - baseline

    fig, ax = plt.subplots(figsize=(14, 4))

    # Normalize colors symmetrically
    max_abs = max(abs(anomaly.min().item()), abs(anomaly.max().item()))
    norm = mcolors.Normalize(vmin=-max_abs, vmax=max_abs)

    colors = plt.cm.get_cmap(cmap)
    years = anomaly.time.dt.year.values

    for i, (year, value) in enumerate(zip(years, anomaly.values)):
        ax.bar(i, 1, width=1, color=colors(norm(value)), edgecolor="none")

    ax.set_xlim(-0.5, len(years) - 0.5)
    ax.set_ylim(0, 1)
    ax.set_xticks([0, len(years)//4, len(years)//2, 3*len(years)//4, len(years)-1])
    ax.set_xticklabels([years[0], years[len(years)//4], years[len(years)//2],
                         years[3*len(years)//4], years[-1]])
    ax.set_yticks([])
    ax.set_title("Temperature Anomaly Stripes", fontsize=14)

    plt.tight_layout()
    return fig

Interactive dashboards with Panel

For exploring multi-scenario, multi-variable output interactively:

import panel as pn
import hvplot.xarray

pn.extension()

class ClimateExplorer:
    """Interactive climate model output explorer."""

    def __init__(self, dataset: xr.Dataset):
        self.ds = dataset
        self.variable = pn.widgets.Select(
            name="Variable",
            options=list(dataset.data_vars),
            value="tas",
        )
        self.scenario = pn.widgets.Select(
            name="Scenario",
            options=["ssp126", "ssp245", "ssp370", "ssp585"],
            value="ssp245",
        )
        self.time_slider = pn.widgets.IntSlider(
            name="Year", start=2020, end=2100, value=2050, step=10,
        )

    @pn.depends("variable.value", "scenario.value", "time_slider.value")
    def map_view(self):
        var = self.variable.value
        decade_start = self.time_slider.value
        data = self.ds[var].sel(
            scenario=self.scenario.value,
            time=slice(f"{decade_start}", f"{decade_start+9}"),
        ).mean("time")

        return data.hvplot.quadmesh(
            x="lon", y="lat",
            projection=ccrs.Robinson(),
            coastline=True,
            cmap="cmocean.thermal" if "tas" in var else "cmocean.rain",
            width=800, height=400,
            title=f"{var} — {self.scenario.value} — {decade_start}s",
        )

    def panel(self):
        return pn.Column(
            pn.Row(self.variable, self.scenario, self.time_slider),
            self.map_view,
        )

# explorer = ClimateExplorer(multi_scenario_ds)
# explorer.panel().servable()

Animated visualizations

from matplotlib.animation import FuncAnimation

def animate_warming(
    annual_data: xr.DataArray,  # (time, lat, lon)
    output_path: str = "warming_animation.mp4",
    fps: int = 5,
):
    """Create an animation of warming over time."""
    fig, ax = plt.subplots(
        figsize=(12, 6),
        subplot_kw={"projection": ccrs.Robinson()},
    )

    # Fixed color scale for consistency across frames
    vmin, vmax = -2, 6
    norm = mcolors.TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)

    def update(frame_idx):
        ax.clear()
        year_data = annual_data.isel(time=frame_idx)
        year = int(annual_data.time.dt.year[frame_idx])

        year_data.plot.pcolormesh(
            ax=ax, transform=ccrs.PlateCarree(),
            cmap=cmocean.cm.balance, norm=norm,
            add_colorbar=False,
        )
        ax.coastlines(linewidth=0.5)
        ax.set_global()
        ax.set_title(f"Temperature Anomaly — {year}", fontsize=14)
        return ax

    anim = FuncAnimation(fig, update, frames=len(annual_data.time), interval=200)
    anim.save(output_path, writer="ffmpeg", fps=fps, dpi=150)
    plt.close()

Tradeoffs

Decision	Option A	Option B
Output format	Static PNG/PDF (publication, reproducible)	Interactive HTML (exploration, web)
Projection	Robinson (balanced world view)	Regional (Lambert/Mercator for specific areas)
Colormap	Sequential (absolute values)	Diverging (anomalies, better for change maps)
Resolution	Native model grid (authentic)	Regridded common grid (comparable across models)
Uncertainty display	Spread shading (intuitive)	Stippling significance (statistically rigorous)

Real-world tools and platforms

• IPCC WGI Interactive Atlas — Built on processed CMIP6 data, allows region/scenario exploration.
• Climate Explorer (KNMI) — Web interface backed by Python processing pipelines.
• Pangeo — Cloud-native climate analysis platform using xarray, Dask, and Jupyter.
• NCL to Python migration — NCAR is actively migrating from NCL (NCAR Command Language) to Python’s GeoCAT library for climate diagnostics.

One thing to remember: Climate visualization is where science meets communication — getting the projection, colormap, and statistical representation right determines whether the graphic accurately informs or inadvertently misleads its audience.