⚙️ How-To Gridded & climate data Joins & relationships Raster & rasterio Gridded & Climate Data · step 11 of 20 Updated

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  47.40 s   190 (52 raised NoDataInBounds)
cell-centre mask, one pass            8.11 s    187
coverage-fraction weights, one pass   3.05 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. The weighted version was fifteen times faster and gave a value to 44 more countries.

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"])
>>> table[["United Kingdom", "Norway", "Singapore"]].sum().round(0)
United Kingdom     961.0
Norway             876.0
Singapore         3614.0

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.

Bar chart of how many of 242 countries contain no grid cell centre at four grid resolutions.
A centre-based method drops these polygons silently or with an exception — and the list grows fast as the grid coarsens.

Step-by-step solution

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

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:

grid countries with no cell centre
0.25° 42
0.5° 52
1° 70
2.5° 101

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 matters. When the real day-one field was written to GeoTIFF instead, the CPC missing value -9.96921e+36 was not declared as nodata, and the mean operation returned −3.23 × 10³¹ mm for Yemen.

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:

country unweighted vs area-weighted
Norway −12.7%
Canada −9.3%
Chile +10.3%
United Kingdom −0.1%

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: 112 of the 242 countries came out as NaN on every day.
  • 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]

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

Bar chart comparing time taken and polygons returned by a clip loop, a centre mask and coverage weights.
The clip loop repeats the geometry work 242 times; the weight matrix does it once and reuses it for all 366 days.

Code examples

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
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
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)
55 polygons have no valid cells at any time: ['Vatican', 'Federated States of Micronesia', ...] ...
no cell centre: 55 polygons; |difference| > 5.0%: 14
Cabo Verde        72.5
Western Sahara    20.6
Puerto Rico       19.4
Palestine         10.9

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.

Explanation

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 47.40 s becomes 3.05 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 12.7% 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.

Two map panels: a small polygon between four cell centres counted by neither, then counted by four partial cells.
The same polygon on the same grid: no centre inside it, but four cells it partly covers.

Edge cases or notes

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

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 3.05 s.

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 12.7% low and Chile's 10.3% 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 112 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.