⚙️ How-To Gridded & climate data Spatial statistics Gridded & Climate Data · step 9 of 20 Updated

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.

Flow from a monthly series through baseline selection and a groupby-month mean to twelve climatology maps, then groupby subtraction to anomalies.
The subtraction matches each month to its own climatology map, so January is compared only with Januaries.

Step-by-step solution

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()
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:

Location Anomaly Baseline s.d. Standardised
Arctic, 80° N 0° E +2.39 °C 4.27 °C +0.56
Equatorial Pacific, 0° 210° E +2.32 °C 1.30 °C +1.78
Indonesia, 0° 115° E +0.99 °C 0.31 °C +3.22

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.

Code examples

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)
monthly_anomalies(air, start="1991-03")
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.

Bar chart of the 2024 global mean temperature anomaly against four baselines, from +0.634 °C for 1991–2020 to +1.095 °C for 1961–1990.
The pattern of anomalies is identical under every baseline; only a constant per cell and month changes.

Explanation

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?".

Two panels contrasting anomalies computed against one long-term mean with anomalies computed against a monthly climatology, with measured values.
A flat baseline turns every winter into a cold anomaly and every summer into a warm one.

Edge cases or notes

  • 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.

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.