How to Extract a Time Series per Polygon from a NetCDF Grid

Problem statement

You have a gridded time series — daily precipitation, monthly temperature, an hourly reanalysis field — and a set of polygons: countries, districts, catchments, sales territories. You want one table with a row per time step and a column per polygon, each value the mean of the grid over that polygon.

The obvious approach is a loop: clip the grid to each polygon, average, repeat. Measured on real data it is both slow and fragile. With the CPC Global Unified Precipitation grid (0.5°, 366 daily fields of 360 × 720 cells for 2024) and the 242 Natural Earth 50 m countries:

method                                time      polygons with a value
rio.clip loop, one polygon at a time  45.05 s   190 (52 raised NoDataInBounds)
cell-centre mask, one pass             0.72 s   187
coverage-fraction weights, one pass    1.21 s   231

The loop stopped on 52 polygons with NoDataInBounds: No data found in bounds. — every one of them a country too small to contain the centre of a single half-degree cell. Both one-pass methods took about a second (Example 4), and the weighted one gave a value to 44 more countries than the centre mask.

Quick answer

Compute, once, what fraction of each grid cell falls inside each polygon. Then every time step is a single sparse matrix product:

import numpy as np
import pandas as pd
import xarray as xr
import geopandas as gpd
import rioxarray  # noqa: F401  (adds the .rio accessor)
import scipy.sparse as sp
from exactextract import exact_extract

pr = xr.open_dataset("precip.2024.nc").precip
pr = pr.assign_coords(lon=((pr.lon + 180) % 360) - 180).sortby("lon")
countries = gpd.read_file("ne_50m_admin_0_countries.zip")[["ADMIN", "geometry"]]

template = xr.ones_like(pr.isel(time=0)).rename(lon="x", lat="y").rio.write_crs(4326)
template.rio.to_raster("template.tif")
cov = exact_extract("template.tif", countries, ["cell_id", "coverage"], output="pandas")

area = np.repeat(np.cos(np.deg2rad(pr.lat.values)), pr.sizes["lon"])
rows = np.concatenate([np.full(len(c), i) for i, c in enumerate(cov["cell_id"])])
cols = np.concatenate(cov["cell_id"]).astype(int)
W = sp.csr_matrix((np.concatenate(cov["coverage"]) * area[cols], (rows, cols)),
                  shape=(len(countries), area.size))

values = pr.values.reshape(pr.sizes["time"], -1)
valid = ~np.isnan(values)
with np.errstate(invalid="ignore"):
    series = (W @ np.where(valid, values, 0).T) / (W @ valid.T)
table = pd.DataFrame(series.T, index=pr.time.values, columns=countries["ADMIN"])
print(table[["United Kingdom", "Norway", "Singapore"]].sum().round(0))

ADMIN
United Kingdom     961.0
Norway             876.0
Singapore         3614.0
dtype: float64

Singapore is the tell: it contains no cell centre at 0.5°, so a clip or a centre mask gives it nothing. Here it gets a value from the one cell it partly covers.

One convention, a coverage matrix, weighted cells

1. Put the grid and the polygons in the same longitude convention

CPC, like most NOAA products, stores longitude from 0 to 360. Natural Earth, like most vector data, uses −180 to 180. Nothing raises an error if they disagree; the Americas simply fall off the grid.

pr = pr.assign_coords(lon=((pr.lon + 180) % 360) - 180).sortby("lon")
print(float(pr.lon[0]), float(pr.lon[-1]))     # -179.75 179.75

-179.75 179.75

Keep latitude descending (north first) as CPC already is. The template raster written in step 3 is then north-up, which is the orientation the cell ids assume.

2. Decide what "inside the polygon" means

There are two defensible definitions, and they disagree for small polygons:

• Cell centre inside the polygon — a cell counts fully or not at all. Simple, and what rio.clip and rasterize do by default.
• Fraction of the cell inside the polygon — every cell the polygon touches counts, weighted by the share it covers.

The centre rule fails on any polygon smaller than a cell. Measured with the same 242 countries:

On a 2.5° reanalysis grid, 42% of the world's countries and territories vanish from a centre-based table.

3. Compute coverage fractions once

exact_extract with the cell_id and coverage operations returns, per polygon, the index of every cell it touches and the fraction of that cell inside it. It only needs the grid geometry, so write a template of ones:

template = xr.ones_like(pr.isel(time=0)).rename(lon="x", lat="y").rio.write_crs(4326)
template.rio.to_raster("template.tif")
cov = exact_extract("template.tif", countries, ["cell_id", "coverage"], output="pandas")

A template of ones also removes a dependency. Written from the real day-one field instead, the GeoTIFF is only safe if the missing value comes with it: here rioxarray carried CPC's missing_value of −9.96921 × 10³⁶ into the file as its nodata, and Yemen's day-one mean came out as a plausible 0 mm (Example 4). A field whose missing value is not declared would feed that number into every mean.

4. Weight each cell by its area

A degree-based cell shrinks towards the poles in proportion to the cosine of its latitude. Leaving that out biases every polygon that spans a range of latitudes. Measured on the 2024 totals:

Norway's error comes from treating its small Arctic cells as if they were as large as its southern ones.

5. Build the matrix and handle missing cells

Stack the weights into a sparse matrix with one row per polygon and one column per cell. Then divide by the weight of the valid cells at each time step, not by the total weight:

valid = ~np.isnan(values)
series = (W @ np.where(valid, values, 0).T) / (W @ valid.T)

CPC is a land-only product with missing values over the sea, so coastal polygons always touch some missing cells. Two common shortcuts both fail, measured:

• Let NaN propagate: 122 of the 242 countries came out as NaN for the year.
• Treat NaN as zero: 12 countries were biased low by more than 10%, the worst South Georgia at −72.9%.

6. Turn the result into a table and aggregate

The product is a polygons × time array; transpose it into a DataFrame indexed by time, and pandas handles the rest:

uk = table["United Kingdom"]
print(uk.resample("MS").sum().round(0).tolist())
# [90.0, 107.0, 87.0, 83.0, 78.0, 40.0, 71.0, 78.0, 94.0, 85.0, 61.0, 87.0]

[90.0, 107.0, 87.0, 83.0, 78.0, 40.0, 71.0, 78.0, 94.0, 85.0, 61.0, 87.0]

The values are area-mean daily depths in millimetres; summing days gives a depth, not a volume.

7. Validate against a centre mask

For large polygons both definitions should agree closely, which is a cheap test that the cell ids and orientation are right. Measured: Russia differed by 0.13%, the United Kingdom by −0.38%. The median absolute difference over all countries was 0.85%; the largest was Cabo Verde at 72.5%, where a single cell centre (93.4 mm) stood in for 16 partly covered cells (54.1 mm).

Daily precipitation per country from CPC

The data used on this page

NOAA's CPC global unified gauge-based daily precipitation for 2024 (0.5°, land only) from NOAA PSL, saved as precip.2024.nc, and Natural Earth's 1:50m countries. PSL drops long transfers, so the setup resumes and retries.

import os
import subprocess

import requests

URLS = ["https://downloads.psl.noaa.gov/Datasets/cpc_global_precip/precip.2024.nc",
        "https://psl.noaa.gov/thredds/fileServer/Datasets/cpc_global_precip/precip.2024.nc"]
for url in URLS * 5:
    if os.path.exists("precip.2024.nc"):
        break
    if subprocess.run(["curl", "-sfL", "-C", "-", "--retry", "3", "-o", "precip.part", url]).returncode == 0:
        os.replace("precip.part", "precip.2024.nc")
if not os.path.exists("ne_50m_admin_0_countries.zip"):
    open("ne_50m_admin_0_countries.zip", "wb").write(requests.get(
        "https://naciscdn.org/naturalearth/50m/cultural/ne_50m_admin_0_countries.zip", timeout=120).content)
print(f"precip.2024.nc {os.path.getsize('precip.2024.nc') / 1e6:.1f} MB, "
      f"countries zip {os.path.getsize('ne_50m_admin_0_countries.zip') / 1e6:.1f} MB")

precip.2024.nc 66.1 MB, countries zip 0.8 MB

Example 1 — building the weight matrix

import numpy as np
import scipy.sparse as sp
import xarray as xr
import rioxarray  # noqa: F401
from exactextract import exact_extract


def polygon_weights(da, polygons, x="lon", y="lat", template_path="template.tif"):
    """Sparse (polygons x cells) matrix of coverage fraction times relative cell area."""
    if float(da[y][0]) < float(da[y][-1]):
        raise ValueError("sort latitude descending first, so the template is north-up")
    grid = xr.ones_like(da.isel({d: 0 for d in da.dims if d not in (x, y)}))
    grid = grid.rename({x: "x", y: "y"}).rio.write_crs(polygons.crs)
    grid.rio.to_raster(template_path)

    cov = exact_extract(template_path, polygons, ["cell_id", "coverage"], output="pandas")
    cell_area = np.repeat(np.cos(np.deg2rad(da[y].values)), da.sizes[x])
    rows = np.concatenate([np.full(len(c), i) for i, c in enumerate(cov["cell_id"])])
    cols = np.concatenate(cov["cell_id"]).astype(np.int64)
    weights = np.concatenate(cov["coverage"]) * cell_area[cols]
    W = sp.csr_matrix((weights, (rows, cols)), shape=(len(polygons), cell_area.size))

    cells = np.diff(W.indptr)
    print(f"{len(polygons)} polygons, {W.nnz:,} weighted cells, "
          f"{int((cells == 0).sum())} polygons touch no cell")
    return W


W = polygon_weights(pr, countries)

242 polygons, 98,875 weighted cells, 0 polygons touch no cell

The matrix depends only on the grid and the polygons. Save it with scipy.sparse.save_npz and every later year of the same product reuses it.

Example 2 — applying it a block of time steps at a time

import numpy as np
import pandas as pd


def polygon_series(da, W, names, time_dim="time", block=60):
    """Area-weighted mean per polygon and time step, a block of time steps at a time."""
    out = []
    n_cells = W.shape[1]
    for start in range(0, da.sizes[time_dim], block):
        chunk = da.isel({time_dim: slice(start, start + block)}).values.reshape(-1, n_cells)
        valid = ~np.isnan(chunk)
        with np.errstate(invalid="ignore", divide="ignore"):
            out.append(((W @ np.where(valid, chunk, 0).T) / (W @ valid.T)).T)
    table = pd.DataFrame(np.vstack(out), index=da[time_dim].values, columns=names)
    empty = table.columns[table.isna().all()]
    if len(empty):
        print(f"{len(empty)} polygons have no valid cells at any time: {list(empty[:5])} ...")
    return table


weighted = polygon_series(pr, W, countries["ADMIN"])

11 polygons have no valid cells at any time: ['American Samoa', 'Saint Helena', 'Pitcairn Islands', 'Bermuda', 'Samoa'] ...

Only one block of time steps is in memory at once, so the same function works on a decade of daily data opened lazily. The eleven empty polygons are small islands where the land-only grid has no valid cell at all — a property of the data, reported rather than hidden.

Example 3 — checking the result against a centre mask

import numpy as np
import scipy.sparse as sp
from rasterio import features
from rasterio.transform import from_origin


def compare_with_centres(da, polygons, table, x="lon", y="lat", threshold=5.0):
    """Annual totals from a cell-centre mask against the coverage-weighted table."""
    res = abs(float(da[x][1] - da[x][0]))
    left, top = float(da[x][0]) - res / 2, float(da[y][0]) + res / 2
    labels = features.rasterize(
        ((geom, i + 1) for i, geom in enumerate(polygons.geometry)),
        out_shape=(da.sizes[y], da.sizes[x]), transform=from_origin(left, top, res, res),
        fill=0, dtype="int32").ravel()
    cols = np.nonzero(labels)[0]
    area = np.repeat(np.cos(np.deg2rad(da[y].values)), da.sizes[x])
    C = sp.csr_matrix((area[cols], (labels[cols] - 1, cols)), shape=(len(polygons), area.size))
    centres = polygon_series(da, C, table.columns)

    total_w, total_c = table.sum(min_count=1), centres.sum(min_count=1)
    diff = ((total_c - total_w) / total_w * 100).dropna()
    print(f"no cell centre: {int(total_c.isna().sum())} polygons; "
          f"|difference| > {threshold}%: {int((diff.abs() > threshold).sum())}")
    return diff.abs().sort_values(ascending=False)


gap = compare_with_centres(pr, countries, weighted)
print(gap.head(4).round(1).to_string())
print(f"median |difference| over {len(gap)} countries: {gap.median():.2f}%")

55 polygons have no valid cells at any time: ['Vatican', 'Federated States of Micronesia', 'Marshall Islands', 'Northern Mariana Islands', 'Guam'] ...
no cell centre: 55 polygons; |difference| > 5.0%: 14
ADMIN
Cabo Verde        72.5
Western Sahara    20.6
Puerto Rico       19.4
Palestine         10.9
median |difference| over 187 countries: 0.85%

The 55 are the 52 polygons with no cell centre plus three whose only centres fall on missing sea cells. The fourteen above 5% are the list to look at before publishing a table built from centres.

Example 4 — the other numbers on this page

import gc
import time
import warnings

from rioxarray.exceptions import NoDataInBounds

warnings.simplefilter("ignore")
names = countries["ADMIN"]

# 1. the three methods, timed
grid = pr.rename(lon="x", lat="y").rio.write_crs(4326)
started, got, no_data = time.perf_counter(), 0, 0
for geom in countries.geometry:
    try:
        grid.rio.clip([geom], countries.crs).mean(["x", "y"]).values
        got += 1
    except NoDataInBounds:
        no_data += 1
    gc.collect()        # each clip leaves its arrays in reference cycles; without this the loop reached 11.6 GB
print(f"rio.clip loop       {time.perf_counter() - started:6.2f} s  {got} polygons ({no_data} NoDataInBounds)")
started = time.perf_counter()
res = 0.5
labels = features.rasterize(((g, i + 1) for i, g in enumerate(countries.geometry)), out_shape=(360, 720),
                            transform=from_origin(-180, 90, res, res), fill=0, dtype="int32").ravel()
centre_cols = np.nonzero(labels)[0]
area = np.repeat(np.cos(np.deg2rad(pr.lat.values)), pr.sizes["lon"])
C = sp.csr_matrix((area[centre_cols], (labels[centre_cols] - 1, centre_cols)), shape=(len(countries), area.size))
centre_table = polygon_series(pr, C, names)
print(f"cell-centre mask    {time.perf_counter() - started:6.2f} s  {int(centre_table.notna().any().sum())} polygons")
started = time.perf_counter()
Wt = polygon_weights(pr, countries)
cover_table = polygon_series(pr, Wt, names)
print(f"coverage weights    {time.perf_counter() - started:6.2f} s  {int(cover_table.notna().any().sum())} polygons")

# 2. countries without a single cell centre, by grid spacing
for res in (0.25, 0.5, 1.0, 2.5):
    lab = features.rasterize(((g, i + 1) for i, g in enumerate(countries.geometry)),
                             out_shape=(int(180 / res), int(360 / res)), transform=from_origin(-180, 90, res, res),
                             fill=0, dtype="int32")
    print(f"{res:4}°: {len(countries) - len(np.unique(lab[lab > 0]))} countries with no cell centre")

# 3. the real day-one field as the template, with the missing value left undeclared
raw = xr.open_dataset("precip.2024.nc", mask_and_scale=False).precip.isel(time=0)
raw = raw.assign_coords(lon=((raw.lon + 180) % 360) - 180).sortby("lon").rename(lon="x", lat="y")
raw.rio.write_crs(4326).rio.to_raster("day1.tif")
yemen = countries[countries["ADMIN"] == "Yemen"]
print(f"Yemen, mean of the undeclared day-one field: {exact_extract('day1.tif', yemen, 'mean', output='pandas')['mean'][0]:.3g} mm")

# 4. area weighting, and two ways to get missing cells wrong
rows = np.concatenate([np.full(len(c), i) for i, c in enumerate(cov["cell_id"])])
cols = np.concatenate(cov["cell_id"]).astype(int)
unweighted = sp.csr_matrix((np.concatenate(cov["coverage"]), (rows, cols)), shape=W.shape)
flat = polygon_series(pr, unweighted, names).sum(min_count=1)
annual = weighted.sum(min_count=1)
for name in ("Norway", "Canada", "Chile", "United Kingdom"):
    print(f"  {name:15} unweighted vs area-weighted {(flat[name] / annual[name] - 1) * 100:+.1f}%")
with np.errstate(invalid="ignore"):
    propagate = pd.Series((W @ values.T).sum(axis=1) / np.asarray(W.sum(axis=1)).ravel(), index=names)
    as_zero = pd.Series((W @ np.where(valid, values, 0).T).sum(axis=1) / np.asarray(W.sum(axis=1)).ravel(), index=names)
print(f"NaN left to propagate: {int(propagate.isna().sum())} countries NaN")
bias = ((as_zero - annual) / annual * 100).dropna()
print(f"NaN treated as zero: {int((bias < -10).sum())} countries more than 10% low; worst {bias.idxmin()} {bias.min():.1f}%")

# 5. Russia, the UK and Cabo Verde by centre and by coverage
for name in ("Russia", "United Kingdom", "Cabo Verde"):
    i = list(names).index(name)
    print(f"  {name:15} centre {centre_table[name].sum():7.1f} mm from {C[i].nnz} cells, "
          f"coverage {annual[name]:7.1f} mm from {Wt[i].nnz} cells, {(centre_table[name].sum() / annual[name] - 1) * 100:+.2f}%")

rio.clip loop        45.05 s  190 polygons (52 NoDataInBounds)
55 polygons have no valid cells at any time: ['Vatican', 'Federated States of Micronesia', 'Marshall Islands', 'Northern Mariana Islands', 'Guam'] ...
cell-centre mask      0.72 s  187 polygons
242 polygons, 98,875 weighted cells, 0 polygons touch no cell
11 polygons have no valid cells at any time: ['American Samoa', 'Saint Helena', 'Pitcairn Islands', 'Bermuda', 'Samoa'] ...
coverage weights      1.21 s  231 polygons
0.25°: 42 countries with no cell centre
 0.5°: 52 countries with no cell centre
 1.0°: 70 countries with no cell centre
 2.5°: 101 countries with no cell centre
Yemen, mean of the undeclared day-one field: 0 mm
11 polygons have no valid cells at any time: ['American Samoa', 'Saint Helena', 'Pitcairn Islands', 'Bermuda', 'Samoa'] ...
  Norway          unweighted vs area-weighted -13.0%
  Canada          unweighted vs area-weighted -9.3%
  Chile           unweighted vs area-weighted +10.1%
  United Kingdom  unweighted vs area-weighted -0.1%
NaN left to propagate: 122 countries NaN
NaN treated as zero: 12 countries more than 10% low; worst South Georgia and the Islands -72.9%
  Russia          centre   464.2 mm from 11722 cells, coverage   463.6 mm from 12934 cells, +0.13%
  United Kingdom  centre   957.2 mm from 133 cells, coverage   960.8 mm from 216 cells, -0.38%
  Cabo Verde      centre    93.4 mm from 1 cells, coverage    54.1 mm from 16 cells, +72.46%

Why a clip loop is slow and cell centres miss

Why the clip loop is slow

Each call to rio.clip rebuilds a mask from the geometry, subsets the array to the polygon's bounds, and then averages all 366 days. The geometry work is repeated for every polygon, and the averaging is repeated inside each iteration.

The weighted method separates the two. Geometry happens once, in exact_extract, and produces 98,875 weights. Time is then a sparse matrix product whose cost is proportional to the number of weights times the number of days. That separation — not a faster clip — is where 45.05 s becomes 1.21 s.

Why cell centres miss small polygons

A cell centre is a point. A polygon smaller than a cell, or a thin one lying between centres, can contain none. The polygon still overlaps the grid; the centre rule simply has nothing to count.

That is why the count of missing countries rises from 42 at 0.25° to 101 at 2.5°. It is also why rio.clip raises NoDataInBounds rather than returning an empty series: after masking there is literally no cell left. Coverage fractions have no such threshold — a polygon covering 3% of one cell gets that cell at 3% weight.

Why the cells must be weighted by area

On a latitude–longitude grid every cell spans the same angles but not the same ground. A cell at 70°N has about a third of the area of one at the equator, since cos 70° ≈ 0.34.

An unweighted mean counts each cell equally, so a country that stretches north over-represents its high-latitude edge. Norway's 13.0% underestimate is exactly that: its wet south coast and its drier Arctic interior count equally per cell although the south covers more ground per cell.

Why missing values have to be renormalised per time step

The denominator of a weighted mean must be the weight of the cells that actually hold data at that time. If a coastal polygon's cells are 30% missing, dividing by the full weight treats the missing 30% as zero rainfall.

On a land-only product the missing cells are fixed, but on satellite or station-derived grids they change from day to day. Dividing by the valid weight at each time step handles both cases with one line.

Overlaps, projected grids and the antimeridian

• Overlapping polygons are fine with weights. A rasterized label grid can hold only one polygon per cell; the sparse matrix gives each overlapping polygon its own row.
• Curvilinear or projected grids need the cell area from the file (areacella in CMIP output) rather than the cosine of latitude.
• Polygons crossing the antimeridian must be split at 180° first; Natural Earth already splits Russia and Fiji.
• Latitude ascending makes rioxarray write a south-up template, and the cell ids no longer line up with a north-first array. Sort descending before building weights.
• Small islands on a land-only grid will stay empty — 11 of 242 here. Report them rather than filling them.
• A centre mask is acceptable when every polygon is many cells wide; check the smallest polygon's cell count before deciding.
• The weights are reusable for any variable on the same grid, and for every year of the same product.
• Units stay the grid's units. Daily precipitation in millimetres averaged over a polygon is still millimetres per day; multiply by area only when you want a volume.

Internal links

• How to clip a NetCDF grid to a polygon in Python — one polygon, and when a clip is enough
• How to take an area-weighted mean over a latitude–longitude grid — the cosine weights on their own
• 0–360 or −180–180: longitude conventions explained — step 1 in depth
• How to calculate zonal statistics in Python — the single-raster version of this task
• How to extract a vegetation index time series for a polygon — the imagery equivalent
• How to rasterize a vector layer in Python — the centre-mask route
• How to resample a gridded time series to monthly or annual values — aggregating the time axis first
• Regridding explained: bilinear, conservative and nearest — coverage fractions are the same idea as conservative regridding
• Raster and vector do not line up in Python — when the polygons land in the wrong cells
• How to extract climate covariates at occurrence points — point sampling of the same climate grids

FAQ

How do I get a time series for each polygon from a NetCDF file?

Compute each polygon's coverage fraction of every grid cell once, multiply by cell area, and apply the resulting sparse matrix to all time steps. For 242 countries and 366 daily fields that took 1.21 s, weights included.

Why does rioxarray clip raise NoDataInBounds for some polygons?

The polygon contains no cell centre, so nothing is left after masking. At 0.5° that happened for 52 of 242 countries; use coverage fractions, or clip with all cells touched, for polygons smaller than a cell.

Do I need to weight grid cells by latitude?

Yes, for any polygon spanning a range of latitudes. Without it Norway's 2024 precipitation came out 13.0% low and Chile's 10.1% high.

How should missing values in the grid be handled?

Divide by the weight of the valid cells at each time step. Letting NaN propagate left 122 of 242 countries empty; treating NaN as zero biased 12 of them by more than 10%.

Is a cell-centre mask ever good enough?

For polygons many cells wide, yes — Russia and the United Kingdom differed by less than 0.4% between the two methods. It breaks down for small or thin polygons, where the difference reached 72.5%.

Can I reuse the weights for another variable or year?

Yes. The matrix depends only on the grid geometry and the polygons, so any field on the same grid can use it without recomputing coverage.