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

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

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 22 ms on the 938 × 73 × 144 array.

Bar chart of the 1991–2020 global mean air temperature: 4.92 °C unweighted against about 14.2 °C with three different area weightings.
Which weighting you choose moves the answer by 0.06 °C. Not weighting moves it by 9.3 °C.

Step-by-step solution

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
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") equals air.mean("lon") everywhere. 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.179 °C

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

7. Sanity-check against the share of area

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.

Code examples

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.

Bar chart of 1979–2025 warming trends: 0.340 °C per decade unweighted, 0.196 weighted, and zonal trends of 1.156 at 80° N, 0.139 at the equator and 0.015 at 60° S.
The unweighted trend is steep because the grid counts the Arctic 80° N row as many times as the equator.

Explanation

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 1979–2025 rise is +1.77 °C. Weighted, the rise is +0.85 °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.

Paired bars comparing the share of grid cells with the share of Earth's area in three latitude bands on a 2.5 degree grid.
Equal-angle grids spend most of their cells where there is least surface.

Edge cases or notes

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

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.

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.