4.7. Vector–Raster Integration

4.7.1. Learning objectives

By the end of this section, you should be able to:

• describe a raster as a grid defined by a CRS, a transform, bands, and a no-data value;

• check that vector and raster inputs are aligned before extracting values;

• extract zonal statistics for polygons and explain which cells each rule includes;

• handle no-data cells, small polygons, and polygons outside the raster;

• choose summaries and resampling methods that fit continuous and categorical data;

• sample a raster along a line; and

• validate extracted values with coverage and sensitivity checks.

Section goal

Combine polygons and lines with raster surfaces without losing track of the grid. A small synthetic raster with known values shows what every extraction rule does before the same steps run on a repository raster.

Vector features describe discrete objects; rasters describe values on a regular grid. Combining them is useful for questions such as:

• What is the mean NDVI [Tucker, 1979] inside each park?

• Which buildings overlap high flood depth?

• What is the elevation along a proposed route?

The answer depends on alignment, resolution, no-data handling, and the unit of analysis. A raster value without its grid definition and CRS is not interpretable.

4.7.2. What defines a raster

A raster is a rectangular grid of cells. Each cell stores one value per band, and the cell represents an area, not a point: a 100 m cell holds one number for 10,000 m². The grid is placed on the Earth by a transform, which gives the coordinates of the top-left corner and the cell size, together with a CRS. Whether a statistic makes sense depends on what the values mean: a continuous raster (elevation, NDVI, temperature) supports means and other numeric summaries, whereas a categorical raster (land cover) stores class codes that must be counted instead. Table 4.26 lists the properties to read before any integration.

Table 4.26 Raster properties and the integration check that each one needs.

Property	What it defines	Integration check
CRS	The coordinate system of the grid	Reproject the vector features to it
Bounds	The extent that the grid covers	Confirm that the features overlap it
Transform and pixel size	Where each cell lies and how large it is	Compare the cell size with the smallest feature
Band data type and meaning	What each value represents and its numeric range	Read the metadata for units, scale factor, and class codes
No-data value	The code that marks missing cells	Mask it, and never average it

The meaning of a value is metadata, not something to guess. Satellite products often store scaled integers, such as a vegetation index multiplied by 10,000, so the metadata decides whether a stored 6,500 means 0.65. The same applies to the no-data code, which is a number like any other until it is masked.

4.7.3. Align the inputs

Before extracting values, compare:

1. CRS;

2. bounds and overlap;

3. pixel size and raster resolution;

4. band data type and meaning; and

5. no-data value and mask.

Reproject the vector features to the raster CRS. Reprojecting vectors only moves vertices, whereas reprojecting a raster resamples its values, so the vector is the layer to move. The cells below build a synthetic raster with known values: a 10 by 10 grid of 100 m cells in a UTM CRS, with values that rise toward the east and a 2 by 2 block of no-data cells standing in for cloud.

Show code cell source

Hide code cell source

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio.io import MemoryFile
from rasterio.mask import mask
from rasterio.transform import from_origin
from shapely.geometry import LineString, box

X0, Y0 = 705500, 5659000        # south-west corner of the toy raster (UTM zone 11N, metres)
NODATA = -9999.0
_keep_alive = []                # in-memory files must outlive the datasets opened from them


def open_raster(array, cell_size, nodata=NODATA, crs='EPSG:32611'):
    # Return an open in-memory GeoTIFF, so the example needs no file on disk.
    rows, cols = array.shape
    memfile = MemoryFile()
    with memfile.open(
        driver='GTiff', height=rows, width=cols, count=1, dtype=array.dtype, crs=crs,
        transform=from_origin(X0, Y0 + rows * cell_size, cell_size, cell_size), nodata=nodata,
    ) as dst:
        dst.write(array, 1)
    _keep_alive.append(memfile)
    return memfile.open()


row_index, col_index = np.mgrid[0:10, 0:10]
ndvi = (0.2 + 0.06 * col_index + 0.01 * row_index).astype('float32')   # rises toward the east
ndvi[2:4, 6:8] = NODATA                                                # a block of no-data cells
src = open_raster(ndvi, 100)


def describe_raster(src):
    return {
        'crs': src.crs.to_string(),
        'bounds': tuple(src.bounds),
        'pixel size (m)': src.res,
        'shape (rows, cols)': (src.height, src.width),
        'bands': src.count,
        'dtype': src.dtypes[0],
        'nodata': src.nodata,
    }


for key, value in describe_raster(src).items():
    print(f'{key:20} {value}')

crs                  EPSG:32611
bounds               (705500.0, 5659000.0, 706500.0, 5660000.0)
pixel size (m)       (100.0, 100.0)
shape (rows, cols)   (10, 10)
bands                1
dtype                float32
nodata               -9999.0

Next come five toy parks that arrive in longitude and latitude, as most vector sources do. Park A covers twelve cell centres, B is smaller than one cell, C overlaps the no-data block, F straddles cell edges, and D lies outside the raster.

Show code cell source

Hide code cell source

def local_box(x0, y0, x1, y1):
    # A box given in metres from the raster's south-west corner.
    return box(X0 + x0, Y0 + y0, X0 + x1, Y0 + y1)


parks = gpd.GeoDataFrame(
    {'park': ['A', 'B', 'C', 'F', 'D']},
    geometry=[
        local_box(110, 610, 490, 890),       # A: twelve cell centres inside
        local_box(110, 110, 140, 140),       # B: smaller than a cell, no centre inside
        local_box(510, 610, 890, 890),       # C: overlaps the no-data block
        local_box(160, 160, 340, 340),       # F: straddles cell edges
        local_box(9000, 9000, 9100, 9100),   # D: outside the raster
    ],
    crs='EPSG:32611',
).to_crs('EPSG:4326')                        # the vector source arrives in longitude and latitude

parks_m = parks.to_crs(src.crs)
print('vector CRS:', parks.crs.to_string(), '-> reprojected to', parks_m.crs.to_string())

overlap = parks_m.geometry.intersects(box(*src.bounds))
print('polygons overlapping the raster:', int(overlap.sum()), 'of', len(parks_m))

vector CRS: EPSG:4326 -> reprojected to EPSG:32611
polygons overlapping the raster: 4 of 5

The same checks on the repository raster [Didan, 2021] follow. The path points to a repository raster; replace it with a local source when working on a different study area. Keep the original file unchanged.

Show code cell source

Hide code cell source

import geopandas as gpd
import rasterio

polygons = gpd.read_file("../data/Parks_Sites_20261019.geojson")

with rasterio.open("../data/MOD13A1_NDVI_20230610.tif") as repo_src:
    if repo_src.crs is None:
        raise ValueError("The raster has no CRS")

    raster_crs = repo_src.crs
    polygons_in_raster_crs = polygons.to_crs(raster_crs)

    print("raster CRS:", repo_src.crs)
    print("raster bounds:", repo_src.bounds)
    print("pixel size:", repo_src.res)

raster CRS: EPSG:3400
raster bounds: BoundingBox(left=131283.8877, bottom=5401931.934, right=864882.3326, top=6693925.173)
pixel size: (470.2554133974359, 526.914045269168)

The raster is in EPSG:3400 with cells of about 470 × 527 m. MOD13A1 is distributed on a 500 m sinusoidal grid, so a grid with non-square cells in a different CRS has been reprojected and resampled. The cell rules below use this stored pixel size, not the product’s nominal 500 m.

4.7.4. Extract values for each polygon

Extracting values for a polygon means choosing which cells belong to it, and Table 4.27 lists the common rules. mask implements the first two. With crop=True it returns only the window around the polygon, and with filled=False it returns a masked array, so cells outside the polygon and no-data cells are masked instead of being replaced by a number. count() on that array is then the number of valid cells.

Table 4.27 Rules for deciding which raster cells belong to a polygon.

Rule	A cell is included when	mask option	Effect
Centre inside	Its centre lies inside the polygon	all_touched=False (default)	A small polygon can include no cell
Any contact	The polygon touches any part of the cell	all_touched=True	Every polygon gets a cell, but edges are over-weighted
Area-weighted	Weighted by the fraction of the cell inside	Not in mask; zonal-statistics tools such as exactextract	Most faithful, and heavier to compute

Two failures follow from these rules. A polygon smaller than a cell may contain no cell centre, so it has zero valid cells, and a polygon that lies completely outside the raster makes mask raise an error. The function below guards against both and reports them instead of stopping.

Show code cell source

Hide code cell source

def summarize(src, polygons, all_touched=False):
    # Zonal mean and valid-cell count for each polygon, keyed by the index.
    if src.closed:
        raise ValueError(
            "The raster dataset is closed. Reopen it before calling summarize()."
        )

    raster_box = box(*src.bounds)
    records = []

    for feature_id, geometry in polygons.geometry.items():
        if not geometry.intersects(raster_box):
            records.append({
                "feature_id": feature_id,
                "valid_cells": 0,
                "mean_value": None,
                "note": "outside raster",
            })
            continue

        values, _ = mask(
            src,
            [geometry],
            crop=True,
            filled=False,
            all_touched=all_touched,
        )

        band = values[0]
        valid = int(band.count())

        records.append({
            "feature_id": feature_id,
            "valid_cells": valid,
            "mean_value": float(band.mean()) if valid else None,
            "note": "" if valid else "no valid cell",
        })

    return pd.DataFrame(records)


# Restore the synthetic raster after the earlier overwrite of src.
src = open_raster(ndvi, 100)

# Match the toy parks to the synthetic raster's CRS.
parks_m = parks.to_crs(src.crs)
by_park = parks_m.set_index("park")

centre_rule = summarize(src, by_park)
print(centre_rule.round(3).to_string(index=False))

feature_id  valid_cells  mean_value           note
         A           12       0.370               
         B            0         NaN  no valid cell
         C            8       0.607               
         F            1       0.390               
         D            0         NaN outside raster

Park A has twelve valid cells. Park C has only eight although it also covers twelve cell centres, because four fall on the no-data block and stay masked. Park B has no valid cell and park D lies outside the raster, and each carries a note that explains its null mean. Changing the rule changes the answer:

Show code cell source

Hide code cell source

any_contact = summarize(src, by_park, all_touched=True)
rules = pd.concat(
    {
        'centre inside': centre_rule.set_index('feature_id')[['valid_cells', 'mean_value']],
        'any contact': any_contact.set_index('feature_id')[['valid_cells', 'mean_value']],
    },
    axis=1,
)
display(rules.round(3))

	centre inside		any contact
	valid_cells	mean_value	valid_cells	mean_value
feature_id
A	12	0.370	12	0.370
B	0	NaN	1	0.340
C	8	0.607	8	0.607
F	1	0.390	9	0.390
D	0	NaN	0	NaN

Parks B and F are the sensitive ones. Under the centre rule B has no cell and F has one, while any contact gives B one cell and F nine. B’s mean goes from missing to 0.340, and F’s mean happens to agree only because the toy values vary smoothly, even though it now rests on nine times as many cells. Parks A and C do not change because their edges follow cell edges. A decision that depends on small or irregular polygons should therefore be tested under both rules.

No-data handling matters as much as the cell rule. Park C shows the cost of treating a no-data code as a measurement:

Show code cell source

Hide code cell source

values, _ = mask(src, [by_park.geometry['C']], crop=True, filled=False)
band = values[0]       # a 3 x 4 window in which every cell belongs to park C

print('masked mean (correct)     :', round(float(band.mean()), 3))
print('no-data replaced by zero  :', round(float(band.filled(0).mean()), 3))
print('no-data kept as -9999     :', round(float(band.data.mean()), 1))

masked mean (correct)     : 0.607
no-data replaced by zero  : 0.405
no-data kept as -9999     : -3332.6

Replacing no-data with zero drags the mean down by a third, and keeping the raw code produces a value that is not an NDVI at all. Neither raises an error.

On the repository raster, the same loop runs over every polygon. This version has no overlap guard, so every polygon must intersect the raster. rasterio.mask.mask clips pixels whose centers or footprints intersect the geometry, depending on the chosen options. The following keeps no-data cells masked:

Show code cell source

Hide code cell source

from rasterio.mask import mask

records = []
with rasterio.open("../data/MOD13A1_NDVI_20230610.tif") as src:
    if raster_crs is None:
        raise ValueError("Run the alignment step before extracting values")
    polygons_in_raster_crs = polygons.to_crs(raster_crs)
    for feature_id, geometry in polygons_in_raster_crs.geometry.items():
        values, _ = mask(src, [geometry], crop=True, filled=False)
        band = values[0]
        valid_cells = int(band.count())
        if valid_cells:
            mean_value = float(band.mean())
            minimum = float(band.min())
            maximum = float(band.max())
        else:
            mean_value = minimum = maximum = None
        records.append(
            {
                "feature_id": feature_id,
                "valid_cells": valid_cells,
                "mean_value": mean_value,
                "minimum": minimum,
                "maximum": maximum,
                "geometry": geometry,
            }
        )

summary = gpd.GeoDataFrame(
    records,
    columns=[
        "feature_id",
        "valid_cells",
        "mean_value",
        "minimum",
        "maximum",
        "geometry",
    ],
    geometry="geometry",
    crs=raster_crs,
)

4.7.5. Categorical and continuous rasters

For a categorical land-cover raster, a mean is usually meaningless. Report class counts or proportions instead. In the toy land-cover raster below, classes 1, 2, and 3 stand for grass, trees, and water, and one polygon spans all three.

Show code cell source

Hide code cell source

landcover = np.ones((10, 10), dtype='int16')
landcover[:, 5:] = 2          # trees in the east
landcover[6:, :] = 3          # water in the south
cover_src = open_raster(landcover, 100, nodata=0)

values, _ = mask(cover_src, [local_box(310, 310, 690, 690)], crop=True, filled=False)
cells_in_polygon = values[0].compressed()

print('mean class (meaningless):', round(float(cells_in_polygon.mean()), 3))
classes, counts = np.unique(cells_in_polygon, return_counts=True)
for cls, count in zip(classes, counts):
    print(f'class {cls}: {count} cells, {count / len(cells_in_polygon):.1%}')

mean class (meaningless): 1.875
class 1: 6 cells, 37.5%
class 2: 6 cells, 37.5%
class 3: 4 cells, 25.0%

The mean of 1.875 describes no land-cover class, whereas the proportions say what is actually there.

4.7.6. Resolution and resampling

Use nearest-neighbor resampling when a new grid must be created from a categorical raster. For a continuous surface, bilinear or cubic resampling can be appropriate, but the choice should be recorded. Table 4.28 lists the common methods.

Table 4.28 Resampling methods and the data they suit.

Method	Best for	Behaviour	Caution
nearest	Categorical rasters	Copies the closest cell value	Blocky, but keeps valid class codes
bilinear	Continuous surfaces	Weighted mean of the four nearest cells	Creates new values and smooths peaks
cubic	Smooth continuous surfaces	Weighted over sixteen cells	Can overshoot the original range
average	Aggregating continuous data to a coarser grid	Mean of the contributing cells	Not for categories
mode	Aggregating categorical data to a coarser grid	Most frequent class	Drops minority classes

Resolution also changes what a polygon sees, even without resampling choices. The cell below aggregates the toy NDVI raster to 200 m cells, ignoring no-data, and extracts the same polygons from both grids.

Show code cell source

Hide code cell source

# Recreate the synthetic 100 m raster.
fine_src = open_raster(ndvi, 100)
by_park = parks.to_crs(fine_src.crs).set_index("park")

# Aggregate 2 × 2 blocks into 200 m cells, ignoring no-data.
blocks = ndvi.reshape(5, 2, 5, 2).swapaxes(1, 2).reshape(5, 5, 4)
is_valid_cell = blocks != NODATA
valid_count = is_valid_cell.sum(axis=2)

coarse = np.where(
    valid_count > 0,
    np.where(is_valid_cell, blocks, 0).sum(axis=2)
    / np.maximum(valid_count, 1),
    NODATA,
).astype("float32")

coarse_src = open_raster(coarse, 200)

# Use fine_src rather than src.
fine = summarize(fine_src, by_park).set_index("feature_id")
coarse_result = summarize(coarse_src, by_park).set_index("feature_id")

comparison = pd.concat(
    {
        "100 m": fine[["valid_cells", "mean_value"]],
        "200 m": coarse_result[["valid_cells", "mean_value"]],
    },
    axis=1,
)

display(comparison.drop(index="D").round(3))

	100 m		200 m
	valid_cells	mean_value	valid_cells	mean_value
feature_id
A	12	0.370	1	0.375
B	0	NaN	0	NaN
C	8	0.607	0	NaN
F	1	0.390	1	0.415

Park A’s mean changes with the grid because each coarse cell averages different neighbours, and park C goes from a mean of 0.607 to no data at all, because its only coarse cell centre lands on the aggregated no-data cell. Neither grid is wrong; they answer at different resolutions, which is why the resolution belongs in the record of the result.

4.7.7. Measurement choices affect the result

Small polygons may contain no complete pixels. Edge cells may be partly covered. Different raster resolutions can produce different summaries for the same polygon. Declare whether you:

• include cells whose centers fall inside the geometry;

• weight edge cells by covered area;

• exclude or impute no-data cells;

• summarize one band or several bands; and

• treat each polygon, pixel, or date as the unit of analysis.

Never replace no-data with zero without checking the source metadata. Zero can be a valid observation, while no-data means that the observation is unavailable.

4.7.8. Quality checks

Besides the structural checks, a coverage ratio shows how well the raster represents each polygon: the area of its valid cells divided by the polygon’s area. A ratio near 1 is reassuring. A ratio of 0 means that nothing could be extracted, a low ratio means that no-data or small size starved the statistic, and a ratio well above 1 means that the polygon is smaller than the cells counted for it.

Show code cell source

Hide code cell source

cell_area = fine_src.res[0] * fine_src.res[1]
qc = centre_rule.assign(polygon_area=parks_m.geometry.area.values)
qc['coverage_ratio'] = qc['valid_cells'] * cell_area / qc['polygon_area']

assert len(centre_rule) == len(parks)
print(qc[['feature_id', 'valid_cells', 'coverage_ratio']].round(2).to_string(index=False))
print('features with zero valid cells:', qc.loc[qc['valid_cells'] == 0, 'feature_id'].tolist())

feature_id  valid_cells  coverage_ratio
         A           12            1.13
         B            0            0.00
         C            8            0.75
         F            1            0.31
         D            0            0.00
features with zero valid cells: ['B', 'D']

Parks B and D have a ratio of 0, F is thin at 0.31, C is limited by no-data at 0.75, and A is above 1 because twelve whole cells cover slightly more than the polygon. The assertions on the repository results are:

Show code cell source

Hide code cell source

assert polygons.crs is not None
assert polygons.geometry.notna().all()
assert len(summary) == len(polygons)
print(summary[["valid_cells", "mean_value"]].describe())

       valid_cells  mean_value
count  3036.000000   282.00000
mean      0.135046  4098.58221
std       0.942259   899.34908
min       0.000000  1571.00000
25%       0.000000  3562.75000
50%       0.000000  4053.50000
75%       0.000000  4558.00000
max      45.000000  6729.00000

Only 282 of the 3,036 parks (9.3%) have any valid cell. Each cell covers about 25 ha (470 × 527 m), which is larger than most parks, so most parks contain no cell centre and their mean is missing, as in toy parks B and D. The means are scaled integers from 1,571 to 6,729, which read as NDVI of 0.16 to 0.67 only if the scale factor is 0.0001, a value to confirm in the product metadata.

Inspect features with zero valid cells, map a sample of extracted values, and compare the result under a second resolution or masking rule when the decision is sensitive to rasterization.

4.7.9. Sampling a raster along a line

A route is a line, so the question “what is the elevation along it?” is answered by sampling the raster at points along the geometry. Sample at an interval no larger than the cell size, or the line can step over whole cells, and treat the no-data code exactly as for polygons. The toy route below runs east across the row that holds the no-data block, sampled every 50 m, half the cell size.

Show code cell source

Hide code cell source

route = LineString([
    (X0, Y0 + 750),
    (X0 + 1000, Y0 + 750),
])

step = 50  # Half the 100 m cell size.
points = [
    route.interpolate(distance)
    for distance in np.arange(0, route.length, step)
]

coordinates = [(point.x, point.y) for point in points]
samples = np.array([
    value[0]
    for value in fine_src.sample(coordinates)
])

is_data = samples != fine_src.nodata

print(
    f"{len(samples)} samples, "
    f"{int((~is_data).sum())} fall on no-data cells"
)
print(
    "mean along the route (no-data masked):",
    round(float(samples[is_data].mean()), 3),
)
print(
    "mean if the no-data code were kept   :",
    round(float(samples.mean()), 1),
)

20 samples, 4 fall on no-data cells
mean along the route (no-data masked): 0.46
mean if the no-data code were kept   : -1999.4

Four samples land on the no-data block, and including them would turn the route’s value into a large negative number. The sampled values can also be summarized as a minimum, maximum, or change between consecutive samples when the question concerns a gradient.

4.7.10. Connection to network analysis

A network edge is a vector feature, so a later network workflow can sample elevation, land cover, or another raster surface along its geometry and assign a documented cost. The same CRS, no-data, and resolution checks apply before using a raster-derived value in Chapter 5 network analysis.

4.7.11. Summary

Vector–raster integration is an alignment and interpretation problem, not merely a file-conversion step. Match CRS and extent, choose the cell rule, resampling, and masking that fit the data type, preserve no-data information, and report the unit and sensitivity of every summary. A coverage ratio, a count of zero-valid-cell features, and a second resolution or masking rule turn those checks into evidence.