Daily precipitation climatology: the annual march of the rain belt¶
GPM IMERG Final (daily) is a global, 0.1° estimate of daily precipitation from NASA’s satellite-and-gauge record. GPM — the Global Precipitation Measurement mission — and IMERG (Integrated Multi-satellitE Retrievals for GPM) merge microwave and infrared observations into one daily rain field; GDEX’s copy spans June 2000 – September 2025.
This notebook reproduces NASA’s IMERG Daily Climatology animation from GDEX data. We
average every calendar day across recent years (2020–2025) into a day-of-year climatology,
smooth it with a 30-day trailing average to remove the noise of isolated heavy-rain days, then
animate the 366-day cycle to watch the Intertropical Convergence Zone (ITCZ) — the equatorial band of heavy tropical rain — and the monsoons march north and south with the seasons.
Compared with the monthly climatology notebook, the daily record resolves seasonal transitions — monsoon onset, the ITCZ’s migration — far more smoothly.
import os
import shutil
import calendar
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.colors import BoundaryNorm
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# Dask
from dask_jobqueue import PBSCluster
from dask.distributed import Client, LocalCluster
import datetime<xarray.core.options.set_options at 0x14e3d197dbd0>Open the dataset¶
IMERG daily is distributed as one NetCDF granule per day. GDEX adds a kerchunk
reference so a single xr.open_dataset(..., engine="kerchunk") call exposes the whole
~25-year daily record as one virtual xarray.Dataset.
The cell below auto-selects the access path: the POSIX reference where /gdex/data is
mounted (NCAR HPC), otherwise the OSDF reference streamed over the network. The OSDF
form’s chunk URLs use the osdf:// scheme, resolved by
pelicanfs.
# POSIX reference where the GDEX data mount is present (NCAR HPC; its chunks live
# under /glade/campaign), otherwise the OSDF-routed reference over the network.
posix_ref = "/gdex/data/d734000/kerchunk/GPM_3IMERGDF_v07.json"
osdf_ref = "https://data.gdex.ucar.edu/d734000/kerchunk/GPM_3IMERGDF_v07-remote-osdf.json"
kerchunk_ref = posix_ref if os.path.exists(posix_ref) else osdf_ref
print(f"Access: {'POSIX mount' if kerchunk_ref == posix_ref else 'OSDF (remote)'}")
print(kerchunk_ref)Access: POSIX mount
/gdex/data/d734000/kerchunk/GPM_3IMERGDF_v07.json
Set up a Dask cluster¶
Building the day-of-year climatology reads several recent years of the daily record, so we farm the work out to a Dask cluster. As in the other GDEX examples, the cluster is chosen automatically by environment:
Casper / Derecho (a PBS scheduler is present) → a
PBSClusterspreads the read across scheduler-managed workers.Laptop / CIRRUS Binder (no PBS) → a
LocalClusteruses the local cores.
def get_pbs_cluster(n_workers=10):
"""PBSCluster on NCAR HPC (Casper/Derecho)."""
scratch = f"/glade/derecho/scratch/{os.environ.get('USER', '')}"
cluster = PBSCluster(
job_name = "imerg-daily",
cores = 1,
memory = "10GiB",
processes = 1,
local_directory = f"{scratch}/dask/spill/",
log_directory = f"{scratch}/dask/logs/",
resource_spec = "select=1:ncpus=1:mem=10GB",
queue = "casper",
walltime = "2:00:00",
interface = "ext",
)
cluster.scale(n_workers)
return cluster
def get_local_cluster():
"""LocalCluster for a laptop or a CIRRUS Binder pod (no PBS scheduler)."""
return LocalCluster()on_hpc = bool(shutil.which("qsub"))
cluster = get_pbs_cluster() if on_hpc else get_local_cluster()
client = Client(cluster)
if on_hpc:
client.wait_for_workers(n_workers=8)
print(f"Detected {'NCAR HPC -> PBSCluster' if on_hpc else 'no PBS -> LocalCluster'}")
clusterBokehUserWarning: reference already known 'p1011'
BokehUserWarning: reference already known 'p1051'
BokehUserWarning: reference already known 'p1289'
BokehUserWarning: reference already known 'p1295'
BokehUserWarning: reference already known 'p1353'
Detected NCAR HPC -> PBSCluster
2026-09-30 12:31:05,722 - distributed.scheduler - ERROR - Couldn't gather keys: {'finalize-hlgfinalizecompute-30cd69f0222f40a2a184bd708ef8e229': 'waiting'}
Day-of-year climatology¶
We average each calendar day (day-of-year 1–366) across recent years (2020–2025) to get the mean annual cycle at daily resolution. Precipitation is already in mm/day, so no unit conversion is needed; we mask the fill value by keeping non-negative rates.
%%time
xr.set_options(use_flox=True)
# Open with LARGE time chunks, not the native 1-day chunks: shrinks the task
# graph ~10x. This matters because dask-jobqueue's scheduler lives in THIS
# kernel, and the native-chunked groupby builds an ~84k-task graph that kills it.
ds = xr.open_dataset(kerchunk_ref, engine="kerchunk", chunks={"time": 120})
precip = ds["precipitation"].where(ds["precipitation"] >= 0) # mm/day; mask fill
precip = precip.sel(time=slice("2020", "2025"))
precip = precip.coarsen(lon=4, lat=4, boundary="trim").mean() # 0.1 -> 0.4 deg
# 'cohorts' is flox's memory-optimal method for day-of-year climatologies.
doy_clim = precip.groupby("time.dayofyear").mean("time", method="cohorts")
doy_clim = doy_clim.transpose("dayofyear", "lat", "lon").compute()
doy_clim.attrs["units"] = "mm/day"
print("shape:", doy_clim.shape)
doy_climshape: (366, 450, 900)
CPU times: user 45.2 s, sys: 4.51 s, total: 49.7 s
Wall time: 2min 32s
Smooth with a 30-day trailing average¶
Individual calendar days are noisy — one big storm in a single year leaves a spot in that day’s mean. Following NASA’s visualization, we apply a 30-day trailing average over the day-of-year axis (wrapping across the December–January boundary) to reveal the smooth seasonal march.
# 30-day trailing mean over day-of-year, wrapping across the year boundary.
n = 30
wrapped = xr.concat(
[doy_clim.isel(dayofyear=slice(-(n - 1), None)), doy_clim], dim="dayofyear")
smooth = (wrapped.rolling(dayofyear=n, center=False).mean().isel(dayofyear=slice(n - 1, None)))
smooth["dayofyear"] = doy_clim["dayofyear"]
smooth.attrs["units"] = "mm/day"
print("smoothed range:", round(float(smooth.min()), 2), "-",
round(float(smooth.max()), 2), "mm/day")
smoothsmoothed range: 0.0 - 50.78 mm/day
Animate the annual cycle¶
Playing the smoothed daily fields as a loop (subsampled to keep the GIF small) reveals the ITCZ and the monsoons marching north and south through the year — the same story as the monthly notebook, but resolved day by day. Following NASA’s visualization, cooler colors are drier and warmer colors wetter (mm/day).
%%time
lat, lon = smooth["lat"], smooth["lon"]
frames = range(0, smooth.sizes["dayofyear"], 3) # every 3rd day keeps the GIF small
fig, ax = plt.subplots(figsize=(8, 4.2),
subplot_kw={"projection": ccrs.PlateCarree()},
constrained_layout=True)
ax.set_global()
ax.coastlines(linewidth=0.4, color="0.5")
mesh = ax.pcolormesh(lon, lat, smooth.isel(dayofyear=0),
transform=ccrs.PlateCarree(), cmap="turbo",
vmin=0, vmax=25, shading="auto", rasterized=True)
cb = plt.colorbar(mesh, ax=ax, orientation="horizontal", shrink=0.7, pad=0.05, extend="max")
cb.set_label("Precipitation (mm/day)")
title = ax.set_title("")
def _label(doy): # day-of-year -> "Mon DD" (2016 is a leap year -> 366 days)
return (datetime.date(2016, 1, 1) + datetime.timedelta(days=int(doy) - 1)).strftime("%b %d")
def update(i):
mesh.set_array(smooth.isel(dayofyear=i).values.ravel())
title.set_text(f"IMERG daily precipitation climatology (30-day rolling mean) - {_label(smooth.dayofyear.values[i])}")
return mesh, title
anim = animation.FuncAnimation(fig, update, frames=frames, interval=120)
anim.save("imerg_daily_climatology.gif", writer=animation.PillowWriter(fps=10), dpi=70)
plt.close(fig)
from IPython.display import Image
Image("imerg_daily_climatology.gif")