⚙️ How-To Raster & rasterio Geometry & topology Raster Analysis · step 9 of 14 Updated

How to Rasterize a Vector Layer in Python

Problem statement

You need a mask: 1 inside the study area, 0 outside, on the same grid as your elevation raster. Or a distance-to-road surface. Or a land-cover grid built from digitised polygons. All of these start by turning vector geometry into cells — and the naive version produces a raster that is the right shape and in the wrong place.

from rasterio.features import rasterize

arr = rasterize([(geom, 1) for geom in gdf.geometry], out_shape=(1000, 1000))

That array is 1,000 × 1,000 of something. It has no transform, so rasterize used the identity transform, and the geometry — with coordinates around 380000, 400000 — fell entirely outside a grid covering 0…1000. The array is all zeros.

Rasterizing correctly means specifying the target grid first and the geometry second. The grid is the decision; the burn is mechanical.

Quick answer 

import rasterio
from rasterio.features import rasterize
import geopandas as gpd

gdf = gpd.read_file("landcover.gpkg")

with rasterio.open("elevation_25m.tif") as template:          # the target grid
    gdf = gdf.to_crs(template.crs)
    burned = rasterize(
        shapes=((geom, value) for geom, value in zip(gdf.geometry, gdf["class_code"])),
        out_shape=(template.height, template.width),
        transform=template.transform,
        fill=0,                       # value for cells no shape covers
        all_touched=False,
        dtype="uint8",
    )
    profile = template.profile.copy()

profile.update(dtype="uint8", count=1, nodata=0, compress="deflate")
with rasterio.open("landcover_25m.tif", "w", **profile) as dst:
    dst.write(burned, 1)
Stack showing the four inputs rasterize needs: shapes with values, output shape, transform and fill value.
Three of the four inputs describe the grid. Only one describes the geometry.
Argument What it decides
out_shape + transform where the grid is and how big its cells are
shapes (geometry, value) pairs — the value is what gets burned
fill the value for cells no geometry covers
all_touched centre-in-polygon (False) or any-touch (True)
dtype must hold your largest burn value — uint8 caps at 255

Step-by-step solution

1. Decide the grid before anything else

There are two ways to get a grid, and choosing the wrong one causes most rasterize problems.

From a template raster — use this whenever the output will be combined with existing raster data:

with rasterio.open("elevation_25m.tif") as tmpl:
    transform, width, height, crs = tmpl.transform, tmpl.width, tmpl.height, tmpl.crs

The result is cell-for-cell aligned with the template, so np.where(mask == 1, elevation, np.nan) is valid. This is almost always what you want.

From the vector's own extent — use this only when the output is standalone:

from rasterio.transform import from_origin

RES = 25
minx, miny, maxx, maxy = gdf.total_bounds
width  = int(np.ceil((maxx - minx) / RES))
height = int(np.ceil((maxy - miny) / RES))
transform = from_origin(minx, maxy, RES, RES)      # note: maxy, and positive RES

from_origin takes the top-left corner — minx, maxy — and positive resolutions, and produces the negative y-scale internally. Passing miny gives a raster that is upside down and 100 km south of the data.

2. Get the CRS right, and make it projected 

gdf = gdf.to_crs(crs)
assert gdf.crs.is_projected, "rasterize in a projected CRS or cell sizes are meaningless"

Rasterizing in EPSG:4326 gives cells measured in degrees, which are neither square nor constant in area. A 0.00025° cell is about 28 m across at the equator and about 17 m at 55° N. Every area, count and distance derived from that grid is wrong by a latitude-dependent factor. See how to choose the right projected CRS.

3. Pick the burn values deliberately

shapes is an iterable of (geometry, value) pairs. Three common patterns:

# (a) a plain mask
shapes = ((geom, 1) for geom in gdf.geometry)

# (b) a coded raster from an attribute
shapes = zip(gdf.geometry, gdf["class_code"].astype("uint8"))

# (c) a value derived per feature
shapes = ((row.geometry, int(row.population / row.area_ha)) for row in gdf.itertuples())

Two rules. The dtype must hold the largest value — burning code 1000 into a uint8 raster wraps it to 232 with no warning. And fill must not collide with a real value: fill=0 is fine when no class is 0, and a disaster when class 0 means "water".

codes = gdf["class_code"]
print(codes.min(), codes.max(), codes.dtype)      # check before choosing dtype

4. Understand what happens where shapes overlap

Two panels showing overlapping polygons burned in sequence, with the later shape overwriting the earlier one.
Last shape wins. Sort deliberately, or the answer depends on row order.

rasterize burns shapes in order, and each one overwrites what is already there. Where two polygons overlap, the value is whichever came later in the iterable.

This makes row order load-bearing, which is a bad property to have by accident:

# make priority explicit: burn low priority first, high priority last
priority = {"water": 3, "urban": 2, "arable": 1, "grassland": 0}
gdf["_prio"] = gdf["class"].map(priority)
gdf = gdf.sort_values("_prio")            # highest priority burned last → wins

If overlaps should not exist at all, that is a data-quality problem to fix before rasterizing — see how to find and fix gaps and overlaps in a polygon coverage.

5. Choose the cell-inclusion rule 

rasterize(shapes, out_shape=..., transform=..., all_touched=False)
  • all_touched=False (default): a cell is burned if its centre is inside the polygon. Cells partition cleanly between neighbouring polygons, so a coverage rasterizes without double-burning.
  • all_touched=True: every cell the geometry touches is burned. Necessary for lines, which otherwise disappear wherever they pass between cell centres.

For lines, all_touched=True is effectively mandatory:

roads = gpd.read_file("roads.gpkg").to_crs(crs)
road_mask = rasterize(((g, 1) for g in roads.geometry),
                      out_shape=(height, width), transform=transform,
                      fill=0, all_touched=True, dtype="uint8")
print(road_mask.sum())      # 48,201 cells   (False gives 11,384 — a broken network)

Points burn to a single cell each regardless of the setting.

Code examples

Example 1: a study-area mask aligned to an existing raster

The workhorse use, in full:

import rasterio, numpy as np, geopandas as gpd
from rasterio.features import rasterize

def mask_from_vector(vector_path, template_path, dst_path, all_touched=False):
    gdf = gpd.read_file(vector_path)
    with rasterio.open(template_path) as tmpl:
        gdf = gdf.to_crs(tmpl.crs)
        arr = rasterize(
            ((geom, 1) for geom in gdf.geometry if geom is not None and not geom.is_empty),
            out_shape=(tmpl.height, tmpl.width),
            transform=tmpl.transform,
            fill=0, all_touched=all_touched, dtype="uint8",
        )
        profile = tmpl.profile.copy()

    profile.update(dtype="uint8", count=1, nodata=None, compress="deflate", tiled=True)
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(arr, 1)
        dst.set_band_description(1, "inside_study_area")

    print(f"{arr.sum():,} of {arr.size:,} cells inside ({100*arr.mean():.1f}%)")
    return arr

mask = mask_from_vector("catchment.gpkg", "elevation_25m.tif", "catchment_mask.tif")
198,447 of 3,038,208 cells inside (6.5%)

The percentage is the check. A mask that is 0.0% or 100.0% means the CRS or the transform is wrong, and it takes one line to notice rather than one afternoon.

Filtering geom is not None and not geom.is_empty prevents rasterize raising on a null geometry — a routine occurrence in real data, and one worth handling at the source instead. See how to remove null and empty geometries.

Example 2: a density surface, and why counts need care

Rasterizing counts is not as simple as burning a value, because burning overwrites rather than accumulates:

import numpy as np, rasterio
from rasterio.features import rasterize

points = gpd.read_file("incidents.gpkg").to_crs(27700)

# WRONG: each point burns 1, so a cell with 20 points still holds 1
naive = rasterize(((g, 1) for g in points.geometry),
                  out_shape=(height, width), transform=transform, dtype="uint16")

# RIGHT: accumulate into the array with the transform's inverse
rows, cols = rasterio.transform.rowcol(transform, points.geometry.x, points.geometry.y)
rows, cols = np.asarray(rows), np.asarray(cols)
inside = (rows >= 0) & (rows < height) & (cols >= 0) & (cols < width)
counts = np.zeros((height, width), dtype="uint16")
np.add.at(counts, (rows[inside], cols[inside]), 1)

print(naive.sum(), counts.sum(), len(points))    # 8412 19003 19003

naive.sum() is the number of occupied cells; counts.sum() is the number of points. The gap is every cell that held more than one. np.add.at accumulates in place and is the standard idiom for this; counts[rows, cols] += 1 does not work, because fancy indexing with duplicate indices applies the update once.

The inside filter matters too: rowcol happily returns negative or out-of-range indices for points outside the grid, and negative indices wrap around to the far edge of the array.

Example 3: rasterizing a folder of layers onto one common grid 

from pathlib import Path
import numpy as np, rasterio, geopandas as gpd
from rasterio.features import rasterize

def build_stack(vector_dir, template_path, dst_path):
    """One band per input layer, all on the template's grid."""
    paths = sorted(Path(vector_dir).glob("*.gpkg"))
    with rasterio.open(template_path) as tmpl:
        shape, transform, crs = (tmpl.height, tmpl.width), tmpl.transform, tmpl.crs
        profile = tmpl.profile.copy()

    bands, names = [], []
    for path in paths:
        gdf = gpd.read_file(path)
        if gdf.empty:
            print(f"  · {path.name:<24} empty, skipped")
            continue
        gdf = gdf.to_crs(crs)
        is_line = gdf.geom_type.str.contains("Line").any()
        arr = rasterize(
            ((g, 1) for g in gdf.geometry if g is not None and not g.is_empty),
            out_shape=shape, transform=transform, fill=0,
            all_touched=is_line,            # lines need it, polygons do not
            dtype="uint8",
        )
        bands.append(arr)
        names.append(path.stem)
        print(f"  ✓ {path.name:<24} {arr.sum():>9,} cells"
              f"{'  (all_touched)' if is_line else ''}")

    profile.update(count=len(bands), dtype="uint8", nodata=None,
                   compress="deflate", tiled=True)
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(np.stack(bands))
        for i, name in enumerate(names, start=1):
            dst.set_band_description(i, name)
    return names

build_stack("layers/", "elevation_25m.tif", "predictors.tif")
  ✓ floodzone.gpkg             412,880 cells
  ✓ railways.gpkg               18,204 cells  (all_touched)
  ✓ roads.gpkg                  48,201 cells  (all_touched)
  · sssi.gpkg                   empty, skipped
  ✓ woodland.gpkg              233,915 cells

Choosing all_touched per layer from its geometry type is the detail that makes this usable unattended: applying one setting to a folder containing both polygons and lines guarantees one of the two is wrong. Band descriptions mean the stack is self-describing later, which matters as much here as in the raster data model.

Explanation

Panels contrasting centre-in-polygon rasterization, which under-includes thin features, with all-touched, which over-includes boundaries.
Two rules, biased in opposite directions. Lines need one; coverages need the other.

Rasterization answers a question with no exact answer: which cells does this shape occupy? A polygon boundary crosses cells, and a cell is atomic — it is in or out, with no partial membership. Every rasterization is therefore a quantisation, and the interesting design decisions are all about how the ambiguity is resolved.

The centre-in-polygon rule (all_touched=False) resolves it by treating each cell as its centre point. This has one excellent property: for a set of non-overlapping polygons, every cell is assigned to at most one of them. A coverage rasterizes to a partition, areas from cell counts are unbiased in aggregate, and the boundary between two polygons is a single line of cells rather than a doubled one. Its cost is under-representation of anything thinner than a cell — which is why it destroys line networks.

The any-touch rule (all_touched=True) resolves it by treating each cell as its full extent. Nothing is missed, so lines stay connected and small polygons survive. But adjacent polygons now share boundary cells, so a coverage no longer partitions, and total area is systematically inflated by roughly the perimeter times half a cell width.

Neither rule is more correct. They are biased in opposite directions, and the right one depends on whether omission or duplication is the worse error for your question.

The second structural point is that rasterization is a lossy, one-way conversion. A polygon has exact coordinates; the raster has a grid of labels. Converting back to vector gives you a stair-stepped approximation of the original, never the original. Round-tripping vector → raster → vector is a coordinate-destroying operation, and the fact that the result looks plausible is what makes it dangerous. Rasterize when the next operation wants a grid — raster algebra, a machine-learning predictor stack, a mask — and keep the vector as the source of record.

Third, the order-dependence of overlapping burns is worth internalising. rasterize is not computing a union or an overlay; it is painting shapes onto a canvas in sequence. Where two shapes overlap, the answer is determined by iteration order, which by default is DataFrame row order, which by default comes from however the file was written. Any correct result from overlapping shapes is an accident unless you sorted deliberately. Sorting by an explicit priority column turns an accident into a decision.

Finally, note the relationship with zonal statistics. Both operations pair polygons with cells; they differ only in direction. Zonal statistics reads cell values grouped by polygon; rasterization writes polygon values into cells. They share the same inclusion rules, and the same choice between partition and cover, for the same reasons.

Edge cases or notes

  • fill must not collide with a real burn value. If class 0 exists, use a different fill and set nodata accordingly.
  • dtype overflow is silent. Burning 300 into uint8 gives 44. Check max() against the dtype.
  • Lines vanish without all_touched=True wherever they pass between cell centres. Always set it for line layers.
  • Points burn to one cell each and do not accumulate. Use np.add.at for counts, as in Example 2.
  • Null or empty geometries raise. Filter them before building the shapes iterable.
  • rasterize takes an iterable, so a generator is fine and avoids materialising millions of tuples.
  • geometry_mask is rasterize with the polarity flipped — it returns a boolean array where True means outside, ready to use as a NumPy mask.
  • rasterio.features.geometry_window limits work to the geometry's window rather than the whole grid, which matters for one small shape on a national grid.
  • gdal_rasterize -a class_code -tr 25 25 -l layer in.gpkg out.tif is the command-line equivalent, and can burn into an existing raster in place.
  • A finer grid then averaging down gives fractional coverage per cell, which is how area weighting in zonal statistics is implemented.

FAQ

Why is my rasterized output all zeros?

No transform was passed, so the identity transform was used and the geometry fell outside the grid. Pass transform= from a template raster or build one with from_origin.

Why did my roads disappear?

Lines only burn cells whose centres they cross. Pass all_touched=True for any line layer.

What happens where polygons overlap?

The shape later in the iterable wins. Sort by an explicit priority column so the outcome is a decision rather than a consequence of row order.

How do I count points per cell?

Not with rasterize — burning overwrites. Convert coordinates to row/col with rasterio.transform.rowcol and accumulate with np.add.at.

What dtype should I use?

The smallest that holds your largest burn value: uint8 up to 255, uint16 up to 65,535. Overflow is silent.

Should the output share a grid with my other rasters?

Yes, whenever you plan to combine them. Take out_shape and transform from a template raster rather than from the vector's bounds.

Can I rasterize into an existing raster instead of a new array?

Yes — pass the existing array as out=. rasterize burns into it in place, leaving cells no shape covers untouched.