How to Calculate a Climatology and Anomalies with xarray

Problem statement

A monthly temperature record is mostly seasonal cycle. In the NCEP/NCAR Reanalysis 1, the area-weighted global mean swings from 12.52 °C in January to 15.98 °C in July every year — a 3.47 °C cycle that has nothing to do with whether a given month was unusual.

To ask "was July 2024 warm?" you need the value expected for July, and the departure from it. The expected value per calendar month is the climatology; the departure is the anomaly. Removing the cycle is not cosmetic: over 1991–2025 the standard deviation of the global monthly series fell from 1.281 °C to 0.310 °C once the climatology was subtracted.

Two shortcuts produce numbers that look reasonable and are wrong. Subtracting one long-term mean instead of a monthly climatology labelled January 2024 at 70° N, 90° E as 17.1 °C below normal; against the proper January climatology it was 0.7 °C below. And grouping daily data by day of year misaligns every date after 28 February in a leap year.

Quick answer

import numpy as np
import xarray as xr

air = xr.open_dataset("air.mon.mean.nc")["air"]

baseline = air.sel(time=slice("1991", "2020"))
climatology = baseline.groupby("time.month").mean("time")      # 12 maps
anomalies = air.groupby("time.month") - climatology            # 938 maps

Measured on the 938-month, 73 × 144 grid: 20 ms for both lines. The 2024 area-weighted global anomaly against 1991–2020 is +0.634 °C in this reanalysis.

[30 30 30 30 30 30 30 30 30 30 30 30]

A complete baseline, a monthly climatology, groupby subtraction

1. Choose a baseline and check that it is complete

A climatology is only as good as its baseline. The World Meteorological Organization's current standard normal is 1991–2020, and a 30-year window is the convention because it averages out year-to-year weather while staying recent.

Count the years that went into each calendar month before using it:

counts = baseline.time.groupby("time.month").count()
print(counts.values)

[30 30 30 30 30 30 30 30 30 30 30 30]

Nothing checks this for you. A baseline that accidentally started at "1991-03" gave January and February 29 years and every other month 30 — and produced a climatology without complaint.

2. Compute the monthly climatology

climatology = baseline.groupby("time.month").mean("time")
print(dict(climatology.sizes))

{'month': 12, 'lat': 73, 'lon': 144}

The time dimension is replaced by month, numbered 1 to 12. For a daily or 6-hourly record the same pattern works with a different key — see step 6.

3. Subtract with groupby arithmetic

anomalies = air.groupby("time.month") - climatology
anomalies = anomalies.drop_vars("month")

xarray pairs every timestep with the climatology map for its calendar month. The result keeps the original time dimension and gains a month coordinate along it, which is worth dropping before concatenating or writing the result — it is redundant with time.

4. Summarise with area weights

A global or regional anomaly series needs cos(latitude) weights, exactly as a mean temperature does:

weights = np.cos(np.deg2rad(anomalies.lat))
global_anomaly = anomalies.weighted(weights).mean(("lat", "lon"))
annual = global_anomaly.groupby("time.year").mean()
for left, right in ((1950, 2016), (1980, 2023), (2000, 2024), (2020, 2025)):
    print(f"{left}  {float(annual.sel(year=left)):+.3f}     {right}  {float(annual.sel(year=right)):+.3f}")

1950  -0.697     2016  +0.485
1980  -0.194     2023  +0.555
2000  -0.240     2024  +0.634
2020  +0.286     2025  +0.508

These are this reanalysis' own values against its own 1991–2020 climate, not an official temperature record. The file also contains two months of 2026, so a naive annual mean for 2026 would be built from January and February alone.

5. Standardise when regions have different variability

A 1 °C anomaly in the tropics is extraordinary; in the Arctic winter it is noise. Dividing by the baseline standard deviation for that calendar month puts every cell on the same footing:

spread = baseline.groupby("time.month").std("time")
z = (anomalies.groupby("time.month") / spread).drop_vars("month")

Measured for January 2024:

Across every month of 2024, 14.9% of cells were more than two standard deviations from normal. In the 1991–2020 baseline itself the share was 4.4%.

6. For daily data, group by calendar day, not day of year

time.dayofyear counts from 1 January, so 1 March is day 61 in a leap year and day 60 otherwise. A day-of-year climatology therefore averages 29 February with 1 March, and day 366 exists only in leap years. Key on the calendar date instead:

files = [f"air.sig995.{year}.nc" for year in range(2020, 2025)]
six_hourly = xr.concat([xr.open_dataset(f)["air"] for f in files], dim="time")
daily = six_hourly.resample(time="D").mean() - 273.15

monthday = daily.time.dt.strftime("%m-%d").rename("monthday")
daily_clim = daily.groupby(monthday).mean()
daily_anom = (daily.groupby(monthday) - daily_clim).drop_vars("monthday")

On 2020–2024 the day-of-year grouping had 5 samples for day 60 but mixed two different calendar dates in it, and only 2 samples for day 366. The two methods disagreed by up to 11.68 °C at a single cell.

7. Re-baseline to compare with a published figure

Anomalies against different baselines differ by a fixed offset per cell and calendar month, so converting between them never needs the raw data again. Measured for 2024:

baseline          2024 global anomaly
1961-1990               +1.095 °C
1981-2010               +0.813 °C
1991-2020               +0.634 °C
full 1948-2025          +0.877 °C

Always state the baseline beside an anomaly. The same year is "+0.63" or "+1.10" depending on that one choice.

Anomalies that refuse an incomplete baseline

The data used on this page

The NCEP/NCAR Reanalysis 1 monthly mean surface air temperature on a 2.5° grid, January 1948 to the latest month, from NOAA PSL (air.mon.mean.nc; the file grows each month, so anything computed up to "now" moves), and the same reanalysis' 6-hourly near-surface air temperature for 2020–2024 (air.sig995.<year>.nc) for the daily examples.

import os

import requests

PSL = "Datasets/ncep.reanalysis/Monthlies/surface/air.mon.mean.nc"
FILES = {
    "oisst-avhrr-v02r01.20240115.nc": ["https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/"
                                       "v2.1/access/avhrr/202401/oisst-avhrr-v02r01.20240115.nc"],
    "air.mon.mean.nc": [f"https://downloads.psl.noaa.gov/{PSL}", f"https://psl.noaa.gov/thredds/fileServer/{PSL}"],
    **{f"air.sig995.{y}.nc": [f"https://downloads.psl.noaa.gov/Datasets/ncep.reanalysis/surface/air.sig995.{y}.nc",
                              f"https://psl.noaa.gov/thredds/fileServer/Datasets/ncep.reanalysis/surface/air.sig995.{y}.nc"]
       for y in range(2020, 2025)},
}
import subprocess

for name, urls in FILES.items():
    for url in urls * 5:                                 # PSL drops connections: resume and retry
        if os.path.exists(name):
            break
        if subprocess.run(["curl", "-sf", "-C", "-", "--retry", "3", "-o", name + ".part", url]).returncode == 0:
            os.replace(name + ".part", name)
    if not os.path.exists(name):
        raise RuntimeError(f"could not download {name}")
print({name: os.path.getsize(name) for name in FILES})

{'oisst-avhrr-v02r01.20240115.nc': 1545085, 'air.mon.mean.nc': 29583703, 'air.sig995.2020.nc': 22123352, 'air.sig995.2021.nc': 22067585, 'air.sig995.2022.nc': 22041061, 'air.sig995.2023.nc': 22084338, 'air.sig995.2024.nc': 22174317}

Example 1 — monthly anomalies that refuse an incomplete baseline

import numpy as np
import xarray as xr


def monthly_anomalies(da, start="1991", end="2020"):
    """Anomalies against a monthly climatology, checking the baseline first."""
    base = da.sel(time=slice(start, end))
    expected = int(end[:4]) - int(start[:4]) + 1
    counts = base.time.groupby("time.month").count()
    short = counts.where(counts < expected, drop=True)
    if short.size:
        detail = ", ".join(f"month {int(m)}: {int(c)} years"
                           for m, c in zip(short.month.values, short.values))
        raise ValueError(f"baseline {start}-{end} is incomplete ({detail})")

    climatology = base.groupby("time.month").mean("time")
    anomalies = (da.groupby("time.month") - climatology).drop_vars("month")
    anomalies.attrs = {**da.attrs, "baseline": f"{start}-{end}",
                       "long_name": f"{da.attrs.get('long_name', da.name)} anomaly"}
    return anomalies, climatology

anomalies, climatology = monthly_anomalies(air)
try:
    monthly_anomalies(air, start="1991-03")
except ValueError as err:
    print(f"ValueError: {err}")

ValueError: baseline 1991-03-2020 is incomplete (month 1: 29 years, month 2: 29 years)

Example 2 — a daily climatology that handles leap years and can be smoothed

def daily_climatology(da, start, end, half_window=0):
    """Calendar-day climatology keyed on month-day, optionally smoothed.

    A 29 February value comes from leap years only. half_window=15 applies
    a 31-day running mean that wraps from 31 December to 1 January.
    """
    base = da.sel(time=slice(start, end))
    key = base.time.dt.strftime("%m-%d").rename("monthday")
    clim = base.groupby(key).mean()
    if half_window:
        n = clim.sizes["monthday"]
        wrapped = xr.concat([clim.isel(monthday=slice(-half_window, None)),
                             clim,
                             clim.isel(monthday=slice(0, half_window))], dim="monthday")
        smoothed = wrapped.rolling(monthday=2 * half_window + 1, center=True).mean()
        clim = smoothed.isel(monthday=slice(half_window, half_window + n))
    samples = base.groupby(key).count().isel({d: 0 for d in da.dims if d != "time"})
    return clim, samples


def daily_anomalies(da, clim):
    key = da.time.dt.strftime("%m-%d").rename("monthday")
    return (da.groupby(key) - clim).drop_vars("monthday")

clim, samples = daily_climatology(daily, "2020", "2024", half_window=15)
print(int(samples.sel(monthday="02-28")), int(samples.sel(monthday="02-29")))

5 2

Smoothing matters more for daily climatologies than monthly ones: each calendar day is an average of only as many values as there are baseline years, so the unsmoothed curve carries day-to-day weather as well as the seasonal cycle.

Example 3 — re-baselining and standardising

def rebaseline(anomalies, start, end):
    """Shift existing anomalies to a new baseline without the raw data."""
    offset = anomalies.sel(time=slice(start, end)).groupby("time.month").mean("time")
    return (anomalies.groupby("time.month") - offset).drop_vars("month")


def standardised_anomalies(da, start="1991", end="2020"):
    """Anomalies divided by the baseline standard deviation of each calendar month."""
    base = da.sel(time=slice(start, end)).groupby("time.month")
    climatology, spread = base.mean("time"), base.std("time")
    anomalies = (da.groupby("time.month") - climatology).drop_vars("month")
    return (anomalies.groupby("time.month") / spread).drop_vars("month")

direct, _ = monthly_anomalies(air, "1961", "1990")
shifted = rebaseline(anomalies, "1961", "1990")
print(float(abs(direct - shifted).max()))

2.002716064453125e-05

Re-baselining the 1991–2020 anomalies to 1961–1990 matched anomalies computed directly from the raw data to 0.00002 °C, which is float32 rounding.

Example 4 — the numbers quoted in the steps

jan = anomalies.sel(time="2024-01").squeeze()
spread = baseline.groupby("time.month").std("time").sel(month=1)
for name, lat, lon in (("Arctic", 80, 0), ("Equatorial Pacific", 0, 210), ("Indonesia", 0, 115)):
    a, sd = float(jan.sel(lat=lat, lon=lon)), float(spread.sel(lat=lat, lon=lon))
    print(f"{name:20} anomaly {a:+.2f}  s.d. {sd:.2f}  standardised {a / sd:+.2f}")
z = standardised_anomalies(air)
print(f"cells beyond 2 s.d.: 2024 {float((abs(z.sel(time='2024')) > 2).mean()):.1%}, "
      f"1991-2020 {float((abs(z.sel(time=slice('1991', '2020'))) > 2).mean()):.1%}")

doy = daily.groupby("time.dayofyear").count().isel(lat=0, lon=0)
print(f"day of year: {int(doy.sel(dayofyear=60))} samples on day 60, {int(doy.sel(dayofyear=366))} on day 366")
by_doy = (daily.groupby("time.dayofyear") - daily.groupby("time.dayofyear").mean()).drop_vars("dayofyear")
print(f"largest disagreement with the calendar-day anomaly: {float(abs(by_doy - daily_anom).max()):.2f} °C")

w = np.cos(np.deg2rad(air.lat))
for start, end in (("1961", "1990"), ("1981", "2010"), ("1991", "2020"), ("1948", "2025")):
    rebased = rebaseline(anomalies, start, end)
    value = float(rebased.sel(time="2024").weighted(w).mean(("lat", "lon")).mean())
    print(f"{start}-{end}   {value:+.3f} °C")

Arctic               anomaly +2.39  s.d. 4.27  standardised +0.56
Equatorial Pacific   anomaly +2.32  s.d. 1.30  standardised +1.78
Indonesia            anomaly +0.99  s.d. 0.31  standardised +3.22
cells beyond 2 s.d.: 2024 14.9%, 1991-2020 4.4%
day of year: 5 samples on day 60, 2 on day 366
largest disagreement with the calendar-day anomaly: 11.68 °C
1961-1990   +1.095 °C
1981-2010   +0.813 °C
1991-2020   +0.634 °C
1948-2025   +0.877 °C

Why the climatology is per calendar month

Why the climatology has to be per calendar month

Subtracting one long-term mean removes the average level but leaves the seasonal cycle in place. Measured against the 1991–2020 all-months mean, the "anomaly" of the global series had a standard deviation of 1.281 °C — identical to the raw series — against 0.310 °C with a monthly climatology.

The distortion is worst where the seasonal cycle is largest. At 70° N, 90° E in Siberia the flat-baseline method called January 2024 17.1 °C below normal. It was 0.7 °C below a normal January. The global mean hides the same effect: January 2024 came out at −1.06 °C and July at +2.31 °C, when the proper anomalies were +0.66 and +0.56.

Why the baseline changes the number but not the pattern

An anomaly against baseline A minus an anomaly against baseline B is the difference between the two climatologies. That difference depends on the cell and the calendar month, never on the year. Measured, the offset between 1961–1990 and 1991–2020 anomalies varied across years by at most 9.5 × 10⁻⁷ °C — floating-point noise.

The offset is not small, and it is not uniform. For July it averaged +0.494 °C globally and ranged from −3.58 °C to +5.69 °C across individual cells. A map of anomalies against an older baseline looks warmer everywhere, and much warmer in some places, without any change in the underlying data.

Why day of year drifts in leap years

dayofyear is a count, not a date. After 28 February, every date in a leap year has a number one higher than the same date in other years, so a day-of-year climatology compares 1 March 2024 with 2 March of the other years and puts 31 December 2024 in a group of its own with 31 December 2020.

On a five-year daily baseline that shift is large, because each calendar day's climatology still contains weather: moving by one day swaps in a different noisy value. Measured, the day-of-year and month-day anomalies for 1 March 2024 differed by up to 6.69 °C at a cell, and for 31 December 2024 by up to 9.84 °C. Non-leap years are affected too — 1 March 2023 differed by 5.59 °C — because their day-60 climatology absorbed two 29 Februaries. A 30-year baseline shrinks the noise; it does not remove the misalignment.

Why standardised anomalies rank places differently

Temperature variability differs by an order of magnitude across the globe. The January standard deviation over 1991–2020 was 4.27 °C in the Arctic at 80° N and 0.31 °C over Indonesia. A +0.99 °C anomaly over Indonesia was therefore 3.22 standard deviations — far rarer than a +2.39 °C anomaly in the Arctic, at 0.56.

Which to map depends on the question. Absolute anomalies answer "how much warmer?"; standardised anomalies answer "how unusual?".

Extra coordinates, silent gaps and partial years

• Groupby arithmetic adds a month coordinate. It is harmless but redundant with time; drop it before concatenating datasets that were grouped differently.
• An incomplete baseline is silent. Starting at "1991-03" left two calendar months with 29 years; check the counts per month.
• Partial final years skew annual summaries. This file ends in February 2026, so an annual mean for 2026 would be a mean of two winter months.
• Mid-month timestamps do not matter here. CMIP6 files stamp months on the 16th; time.month is unaffected.
• Non-standard calendars work. On a noleap cftime axis time.month behaves normally, and day of year never reaches 366.
• Precipitation anomalies are often relative. Dividing by the climatology, rather than subtracting it, compares wet and dry regions more fairly.
• Kelvin and Celsius anomalies are the same number. A difference of 1 K is a difference of 1 °C, so an anomaly needs no unit conversion.
• Large Dask-backed grids benefit from flox. xarray uses it automatically for groupby reductions when it is installed; it was not installed for these measurements.

Internal links

• How to take an area-weighted mean over a latitude–longitude grid — the weights behind every global anomaly number
• How to resample a gridded time series to monthly or annual values — building the monthly series a climatology needs
• How to select a time range and location from an xarray Dataset — slicing the baseline period
• How to open a NetCDF file in Python with xarray — loading the monthly file
• How to open hundreds of NetCDF files as one dataset — daily records spread across yearly files
• Fixing cftime and out-of-bounds datetime errors in xarray — climate-model calendars
• How to extract a time series per polygon from a NetCDF grid — regional anomaly series
• Lazy loading explained — why these groupby operations can run on data larger than memory

FAQ

What is the difference between a climatology and an anomaly?

A climatology is the expected value for each calendar month or day, averaged over a baseline period. An anomaly is an observed value minus the climatology for its calendar month or day.

Which baseline period should I use?

1991–2020 is the current World Meteorological Organization standard normal. Whatever you choose, state it: the 2024 global anomaly here was +0.634 °C against 1991–2020 and +1.095 °C against 1961–1990.

Why does my anomaly map still show a seasonal pattern?

Usually because a single long-term mean was subtracted rather than a monthly climatology. That left the seasonal cycle intact and labelled a normal Siberian January as 17.1 °C below normal.

How do I compute daily anomalies with leap years?

Group by the month-day string from strftime rather than by day of year. Day of year shifts every date after 28 February in leap years and creates a day 366 that only leap years contribute to.

What is a standardised anomaly?

The anomaly divided by the baseline standard deviation for that calendar month. It measures how unusual a value is: a +0.99 °C anomaly over Indonesia was 3.22 standard deviations, a +2.39 °C Arctic anomaly only 0.56.

Can I change the baseline without recomputing everything?

Yes. Subtract the mean of the existing anomalies over the new baseline period, per calendar month. The offset is a constant per cell and month, so no raw data is needed.