How to Take an Area-Weighted Mean over a Latitude–Longitude Grid

On this page

Problem statement

A global mean of a gridded field looks like a one-liner:

air.mean(("lat", "lon"))

On a regular latitude–longitude grid it is wrong, and not by a rounding error. Every row of a 2.5° grid holds 144 cells, whether that row runs around the equator or circles the pole in a few hundred kilometres. A plain mean gives every cell the same vote, so the small polar cells outvote the large tropical ones.

Measured on the NCEP/NCAR Reanalysis 1 monthly surface air temperature, averaged over 1991–2020:

unweighted mean            4.92 °C
cos(latitude) weighted    14.24 °C

That is a 9.3 °C error in a single number. It distorts change over time as well: the unweighted 1979–2025 trend came out at 0.340 °C per decade against 0.196 °C per decade weighted — 73% too steep — because the region warming fastest, the Arctic, is exactly the region the unweighted mean over-counts.

Quick answer

import numpy as np
import xarray as xr

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

weights = np.cos(np.deg2rad(air.lat))                      # a DataArray on "lat"
global_mean = air.weighted(weights).mean(("lat", "lon"))   # one value per month

The weights need only the latitude dimension; xarray broadcasts them across longitude and time. The result keeps time: 938 monthly global means, computed in 21 ms on the 938 × 73 × 144 array.

Which weighting you choose moves the answer by 0.06 °C. Not weighting moves it by 9.3 °C.

What you get — from 1. Confirm the grid is regular in latitude and longitude

[90.  87.5 85. ] [-85.  -87.5 -90. ]
[-2.5] [2.5]

Check the grid, weight by latitude, reduce with NaNs

Only the first column is right; the third is wrong and looks reasonable.

1. Confirm the grid is regular in latitude and longitude

Weighting by cos(latitude) is correct for cells bounded by lines of constant latitude and longitude with uniform spacing. Check before trusting it:

print(air.lat.values[:3], air.lat.values[-3:])
print(np.unique(np.diff(air.lat.values)), np.unique(np.diff(air.lon.values)))

[90.  87.5 85. ] [-85.  -87.5 -90. ]
[-2.5] [2.5]

Descending latitude makes no difference to the weights. A curvilinear model grid with two-dimensional lat and lon arrays, or a rotated-pole regional grid, is a different case: use the cell areas the model publishes rather than this recipe.

2. Build the weights from latitude, in radians

weights = np.cos(np.deg2rad(air.lat))
weights.name = "weights"

The mistake to avoid is forgetting the conversion. np.cos(air.lat) treats 60 as 60 radians and produces negative weights for some rows. Measured, the resulting "global mean" was −24.0 °C; with the absolute value of those weights it was 4.84 °C — close enough to the unweighted answer that nobody would question it.

3. Let `.weighted()` handle missing values

For a field with gaps — sea-surface temperature, where land is NaN, or a grid clipped to a region — the denominator must count only the cells that hold data. .weighted() does that. Measured on one day of NOAA OISST v2.1 at 0.25°, where 33.3% of cells are land:

sst = xr.open_dataset("oisst-avhrr-v02r01.20240115.nc")["sst"].squeeze(drop=True)
w = np.cos(np.deg2rad(sst.lat))

right = sst.weighted(w).mean(("lat", "lon"))
wrong = (sst * w).mean() / w.mean()         # denominator includes land rows

careful = (sst * w).sum() / (w * sst.notnull()).sum()
print(f"land cells {float(sst.isnull().mean()):.1%}")
for name, value in [("unweighted", sst.mean()), (".weighted(cos lat)", right),
                    ("manual, NaN-aware denominator", careful), ("manual, every cell in denominator", wrong)]:
    print(f"{name:36} {float(value):.3f} °C")

land cells 33.3%
unweighted                           14.140 °C
.weighted(cos lat)                   18.818 °C
manual, NaN-aware denominator        18.818 °C
manual, every cell in denominator    20.159 °C

The hand-rolled version was 1.34 °C too warm, because the mean of the weights over every cell is not the mean over ocean cells.

4. Reduce over the horizontal dimensions only

Name the dimensions. mean(("lat", "lon")) returns a time series; mean() with no arguments averages time as well and returns one number.

A zonal mean — averaging along longitude only — needs no weights at all, because every cell in a row has the same area. Verified: air.weighted(w).mean("lon") and air.mean("lon") differ by at most 0.000015 °C anywhere, which is float32 rounding. Weights start to matter at the moment latitudes are combined.

5. Select the region first, then weight

A regional box needs the same treatment, and the error grows with the range of latitudes inside it:

air180 = air.assign_coords(lon=((air.lon + 180) % 360) - 180).sortby("lon")
europe = air180.sel(lat=slice(70, 35), lon=slice(-10, 40))
europe_mean = europe.weighted(np.cos(np.deg2rad(europe.lat))).mean(("lat", "lon"))

Measured on the 1991–2020 mean:

Region	Unweighted	Weighted	Error
Tropics, 20° S–20° N	25.373 °C	25.385 °C	−0.012 °C
Europe, 35–70° N, 10° W–40° E	9.253 °C	10.428 °C	−1.175 °C
Arctic, 60–90° N	−10.077 °C	−7.099 °C	−2.977 °C
Northern Hemisphere	7.817 °C	15.323 °C	−7.506 °C

The latitude slice runs from high to low because this file stores latitude descending — see selecting a time range and location.

6. Decide whether cos(latitude) is precise enough

cos(latitude) is the area of a cell per unit of latitude at its centre. The exact spherical area uses the cell edges, and an ellipsoid adds the flattening of the Earth. Measured on the 1991–2020 global mean:

cos(lat)                  14.239 °C
exact spherical band      14.229 °C
WGS84 ellipsoid           14.177 °C

The three differ by 0.06 °C. For a published global number, use exact band areas; for everything else, cos(latitude) is fine.

The whole problem is visible in one comparison. On this grid the rows at or poleward of 60° hold 35.6% of the cells and 14.5% of the area. Any statistic that treats cells as equal gives that region two and a half times its real influence.

Weighted global and regional temperature means

Example 2 reports the coverage alongside the mean so the denominator is never invisible.

The data used on this page

Two NOAA files. NCEP/NCAR Reanalysis 1 monthly mean surface air temperature on a 2.5° grid (air.mon.mean.nc, January 1948 to the latest month, from NOAA PSL; the monthly file grows each month, so the month count and anything computed up to "now" moves), and one day of NOAA OISST v2.1 at 0.25°, 15 January 2024, from NCEI.

import os

import numpy as np
import requests
import xarray as xr

FILES = {
    "air.mon.mean.nc": "https://psl.noaa.gov/thredds/fileServer/Datasets/ncep.reanalysis/Monthlies/surface/air.mon.mean.nc",
    "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",
}


def download(url, name, attempts=10):
    """PSL's server often drops the connection part-way; resume from the last byte rather than restart."""
    part = name + ".part"
    for _ in range(attempts):
        have = os.path.getsize(part) if os.path.exists(part) else 0
        try:
            with requests.get(url, headers={"Range": f"bytes={have}-"} if have else {}, stream=True, timeout=600) as r:
                if r.status_code == 416:                 # nothing left to fetch
                    break
                r.raise_for_status()
                with open(part, "ab" if r.status_code == 206 else "wb") as out:
                    for chunk in r.iter_content(1 << 20):
                        out.write(chunk)
                total = int(r.headers.get("Content-Range", "/0").split("/")[-1] or 0) or int(r.headers["Content-Length"])
            if os.path.getsize(part) >= total:
                break
        except requests.RequestException:
            continue
    else:
        raise RuntimeError(f"could not download {name}")
    os.replace(part, name)


for name, url in FILES.items():
    if not os.path.exists(name):
        download(url, name)

air = xr.open_dataset("air.mon.mean.nc")["air"].load()
sst = xr.open_dataset("oisst-avhrr-v02r01.20240115.nc")["sst"].squeeze(drop=True).load()
print(air.shape, str(air.time.values[0])[:7], str(air.time.values[-1])[:7], air.attrs.get("units"))

(938, 73, 144) 1948-01 2026-02 degC

Example 1 — weights for any regular latitude–longitude grid

import numpy as np
import xarray as xr
from pyproj import Geod


def cell_area_weights(lat, method="cos"):
    """Relative area of each latitude row of a regular lat-lon grid.

    method="cos"        cos(latitude) at the cell centre
    method="band"       exact spherical area between the cell edges
    method="ellipsoid"  WGS84 band areas from pyproj
    """
    phi = lat.values.astype("float64")
    if method == "cos":
        values = np.cos(np.deg2rad(phi)).clip(min=0)
    else:
        half = abs(phi[1] - phi[0]) / 2
        upper = np.clip(phi + half, -90, 90)
        lower = np.clip(phi - half, -90, 90)
        if method == "band":
            values = np.sin(np.deg2rad(upper)) - np.sin(np.deg2rad(lower))
        elif method == "ellipsoid":
            geod = Geod(ellps="WGS84")
            values = np.array([
                abs(geod.polygon_area_perimeter([0, 1, 1, 0], [lo, lo, hi, hi])[0])
                for lo, hi in zip(lower, upper)
            ])
        else:
            raise ValueError(f"unknown method {method!r}")
    return xr.DataArray(values / values.max(), coords=lat.coords, dims=lat.dims,
                        name=f"area_weight_{method}")

It works in float64 on purpose. With float32 latitudes, np.deg2rad of 90 lands a hair past π/2 and the pole row gets a weight of −4.4 × 10⁻⁸ — harmless in size, but a negative weight is never what you meant. The clip removes it.

Example 2 — a NaN-aware regional or global mean that reports its coverage

def area_weighted_mean(da, weights=None, lat="lat", lon="lon"):
    """Mean over the horizontal dimensions, keeping every other dimension.

    Returns the mean and the fraction of the region's area that had data,
    so a mean built from a handful of cells is visible as such.
    """
    if weights is None:
        weights = cell_area_weights(da[lat])
    mean = da.weighted(weights).mean((lat, lon))
    total = float(weights.sum()) * da.sizes[lon]
    coverage = (da.notnull() * weights).sum((lat, lon)) / total
    return mean, coverage

mean, coverage = area_weighted_mean(sst)
print(f"{float(mean):.3f} °C over {float(coverage):.1%} of the globe's area")

18.818 °C over 71.4% of the globe's area

Sea surface covers about 71% of the Earth, so the coverage figure doubles as a check that the land mask and the weights agree.

Example 3 — a report that shows what weighting changes

def weighting_report(da, start="1979", end="2025", lat="lat", lon="lon"):
    """Mean and linear trend of the annual series, with and without weights."""
    annual = da.sel(time=slice(start, end)).groupby("time.year").mean()
    weights = cell_area_weights(da[lat])
    series = {
        "unweighted": annual.mean((lat, lon)),
        "area-weighted": annual.weighted(weights).mean((lat, lon)),
    }
    years = annual.year.values
    for name, values in series.items():
        slope = np.polyfit(years, values.values, 1)[0] * 10
        print(f"{name:14} mean {float(values.mean()):7.3f}   trend {slope:+.3f} per decade")
    return series

_ = weighting_report(air)

unweighted     mean   4.839   trend +0.340 per decade
area-weighted  mean  14.209   trend +0.196 per decade

Run it on any new gridded dataset before building on its global numbers. If the two rows agree closely, the field has little latitude structure; if they do not, every unweighted statistic downstream is suspect.

The unweighted trend is steep because the grid counts the Arctic 80° N row as many times as the equator.

Example 4 — the other numbers on this page

import time

clim = air.sel(time=slice("1991", "2020")).mean("time")
w = np.cos(np.deg2rad(air.lat))
print(f"1991-2020: unweighted {float(clim.mean()):.2f}, cos(lat) {float(clim.weighted(w).mean()):.2f}")
print(f"forgot deg2rad: {float(clim.weighted(np.cos(air.lat)).mean()):.1f}; "
      f"abs() of those weights: {float(clim.weighted(abs(np.cos(air.lat))).mean()):.2f}")
for method in ("cos", "band", "ellipsoid"):
    print(f"{method:10} {float(clim.weighted(cell_area_weights(air.lat, method)).mean()):.3f}")

start = time.perf_counter()
series = air.weighted(w).mean(("lat", "lon"))
elapsed = time.perf_counter() - start
print(f"{series.size} monthly means from {air.shape}")
print(f"zonal mean, weighted minus unweighted: at most {float(abs(air.weighted(w).mean('lon') - air.mean('lon')).max()):.6f} °C")

c180 = clim.assign_coords(lon=((clim.lon + 180) % 360) - 180).sortby("lon")
for name, box in [("tropics", c180.sel(lat=slice(20, -20))),
                  ("Europe", c180.sel(lat=slice(70, 35), lon=slice(-10, 40))),
                  ("Arctic", c180.sel(lat=slice(90, 60))),
                  ("NH", c180.sel(lat=slice(90, 0)))]:
    u, v = float(box.mean()), float(box.weighted(np.cos(np.deg2rad(box.lat))).mean())
    print(f"{name:8} unweighted {u:8.3f}  weighted {v:8.3f}  error {u - v:+.3f}")

band = cell_area_weights(air.lat, "band")
high = abs(air.lat) >= 60
print(f"rows at or beyond 60°: {float(high.mean()):.1%} of cells, "
      f"{float(band.where(high).sum() / band.sum()):.1%} of area; "
      f"north pole row {float(band.sel(lat=90) / band.sum()):.3%} of the globe")
print(f"float32 pole weight: {np.cos(np.deg2rad(np.float32(90))):.1e}")

annual = air.sel(time=slice("1979", "2025")).groupby("time.year").mean()
years = annual.year.values
for lat in (80, 0, -60):
    print(f"zonal trend at {lat:+d}: {np.polyfit(years, annual.sel(lat=lat).mean('lon'), 1)[0] * 10:.3f} per decade")
for name, s_ in (("unweighted", annual.mean(("lat", "lon"))), ("weighted", annual.weighted(w).mean(("lat", "lon")))):
    print(f"{name} rise 1979-2025 (fitted): {np.polyfit(years, s_, 1)[0] * (years[-1] - years[0]):+.2f}")

print(f"weighted median {float(clim.weighted(w).quantile(0.5)):.1f}, unweighted {float(clim.median()):.1f}")
for bad in (w.values, w.where(w > 0.5)):
    try:
        clim.weighted(bad)
    except Exception as err:
        print(type(err).__name__, "-", str(err).splitlines()[0][:90])

1991-2020: unweighted 4.92, cos(lat) 14.24
forgot deg2rad: -24.0; abs() of those weights: 4.84
cos        14.239
band       14.229
ellipsoid  14.177
938 monthly means from (938, 73, 144)
zonal mean, weighted minus unweighted: at most 0.000015 °C
tropics  unweighted   25.373  weighted   25.385  error -0.012
Europe   unweighted    9.253  weighted   10.428  error -1.175
Arctic   unweighted  -10.077  weighted   -7.099  error -2.977
NH       unweighted    7.817  weighted   15.323  error -7.506
rows at or beyond 60°: 35.6% of cells, 14.5% of area; north pole row 0.012% of the globe
float32 pole weight: -4.4e-08
zonal trend at +80: 1.156 per decade
zonal trend at +0: 0.139 per decade
zonal trend at -60: 0.015 per decade
unweighted rise 1979-2025 (fitted): +1.56
weighted rise 1979-2025 (fitted): +0.90
weighted median 19.2, unweighted 8.8
ValueError - `weights` must be a DataArray
ValueError - `weights` cannot contain missing values. Missing values can be replaced by `weights.fillna

The weighted mean over the whole array, timed on this machine:

print(f"{elapsed * 1000:.0f} ms for {series.size} weighted global means")

21 ms for 938 weighted global means

Why degree cells shrink towards the poles

Why a degree cell shrinks towards the poles

On a sphere of radius R, the area between latitudes φ₁ and φ₂ and longitudes λ₁ and λ₂ is R²(λ₂ − λ₁)(sin φ₂ − sin φ₁). For a narrow cell, sin φ₂ − sin φ₁ is very nearly cos φ multiplied by the cell height, so area is proportional to cos φ.

A 2.5° cell at 60° therefore has half the area of one at the equator, and a cell centred at 87.5° has 4.4% of it. The grid stores the same number of values in each row, so the data density per square kilometre rises steeply towards the poles.

Why the unweighted mean is so much colder

Because the over-represented cells are also the coldest ones. The 35.6% of cells at or beyond 60° cover 14.5% of the planet and include both polar regions, where the 1991–2020 mean temperature is far below zero.

The unweighted average is really an average over rows of latitude. It answers the question "what is the mean temperature of a randomly chosen grid cell", which no one asks.

Why the trend is wrong as well

A weighting error that was constant in time would shift the mean and leave the trend alone. This one does not, because different latitudes warmed at different rates. Measured zonal-mean trends over 1979–2025: 1.156 °C per decade at 80° N, 0.139 at the equator and 0.015 at 60° S.

The unweighted series gives the fast-warming Arctic rows the same weight as the tropics, so its trend is 0.340 °C per decade and its fitted 1979–2025 rise is +1.56 °C. Weighted, the rise is +0.90 °C. Any comparison of warming between datasets or regions has to weight first.

Why the pole rows need a moment's thought

This grid has rows exactly at 90° N and 90° S. Each "cell" there is really a cap of 1.25° around the pole, with a small but non-zero area: 0.012% of the globe for the northern row, using exact band weights. cos(90°) gives zero (or a tiny negative number in float32), so the pole values are ignored.

Here the difference is negligible, which is why cos(latitude) and exact bands agreed to 0.01 °C. On a grid whose rows sit at cell edges rather than centres, the half-cells at each end matter more, and exact band weights are the safer choice.

Why the ellipsoid barely matters

The Earth is flattened by about one part in 298, which moves area slightly away from the poles relative to a sphere. Measured, WGS84 band areas changed the 1991–2020 global mean by 0.05 °C against the spherical band. That is well inside the uncertainty of any reanalysis, and two orders of magnitude smaller than the error from not weighting.

Equal-angle grids spend most of their cells where there is least surface.

Weight types, Gaussian grids and published areas

Weights must be a DataArray. Passing a NumPy array raises ValueError saying the weights must be a DataArray; wrap them with xr.DataArray(values, coords={"lat": lat}).
Weights cannot contain NaN. The error message suggests weights.fillna(0), which is right when the NaN marks a cell that should not count.
Gaussian grids are not evenly spaced in latitude. Spectral models output rows at Gaussian latitudes; use band areas with edges halfway between rows rather than cos(latitude) alone.
Curvilinear and rotated grids need published areas. CMIP6 supplies areacella for atmosphere grids and areacello for ocean grids; weight by those.
Land-only means need a land fraction. Multiply the area weight by the land fraction (sftlf in CMIP6) so coastal cells count in proportion.
Weighted statistics go beyond the mean. .weighted(w).std(), .sum() and .quantile() exist; the weighted median of the 1991–2020 field was 19.2 °C against 8.8 °C unweighted.
Tropical boxes barely change. Over 20° S–20° N the error was 0.012 °C; the correction matters in proportion to the span of latitudes averaged.
Dask-backed arrays work unchanged. The weights are a small in-memory array on one dimension and broadcast lazily.

Internal links

How to calculate a climatology and anomalies with xarray — a global anomaly series needs these weights
How to clip a NetCDF grid to a polygon in Python — regional means where the NaN handling matters
How to extract a time series per polygon from a NetCDF grid — weighted means for many regions at once
0–360 or −180–180: longitude conventions explained — converting the grid before selecting a regional box
How to select a time range and location from an xarray Dataset — slicing a descending latitude axis
Regridding explained — conservative regridding is built on the same cell areas
How to aggregate millions of points into a grid with DuckDB — the same degree-cell distortion, for point counts
Projected vs geographic CRS — why a degree is not a unit of area

FAQ

Why do I need to weight a global mean by latitude?

Because the cells of a latitude–longitude grid shrink towards the poles while each row holds the same number of values. Unweighted, the 1991–2020 global mean air temperature came out at 4.92 °C instead of 14.24 °C.

Should I use cos(latitude) or the exact cell area?

cos(latitude) for almost everything. Exact spherical bands changed the global mean by 0.01 °C and a WGS84 ellipsoid by a further 0.05 °C; skipping weights entirely changed it by 9.3 °C.

Does xarray's weighted mean skip NaN values?

Yes. The denominator counts only cells with data, which is why it matched a careful manual calculation on sea-surface temperature, while a manual version that divided by every cell's weight was 1.34 °C too warm.

Do I need weights for a zonal mean?

No. Every cell along a row of latitude has the same area, so a mean over longitude alone is identical with or without weights. Weights matter once latitudes are combined.

How do I weight a curvilinear or rotated grid?

Use the cell-area variable distributed with the model output, such as areacella in CMIP6, rather than cos(latitude). The latitude of a cell no longer determines its area on those grids.

Does area weighting change trends as well as averages?

Yes, whenever warming varies with latitude. Over 1979–2025 the unweighted trend was 0.340 °C per decade and the weighted one 0.196, because the Arctic warmed at 1.156 °C per decade.