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.

import geopandas as gpd
from rasterio.features import rasterize

gdf = gpd.read_file("leeds.gpkg")                 # the Leeds district, British National Grid
arr = rasterize([(geom, 1) for geom in gdf.geometry], out_shape=(1000, 1000))
print(arr.shape, "cells burned:", int(arr.sum()), "| geometry x from", round(gdf.total_bounds[0]))

(1000, 1000) cells burned: 0 | geometry x from 413216

That array is 1,000 × 1,000 of something. It has no transform, so rasterize used the identity transform, and the geometry — with coordinates in the hundreds of thousands of metres — 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 numpy as np
import rasterio

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

with rasterio.open("elevation_30m.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_30m.tif", "w", **profile) as dst:
    dst.write(burned, 1)
print(f"{len(gdf):,} polygons → {(burned > 0).sum():,} of {burned.size:,} cells burned "
      f"with {len(np.unique(burned)) - 1} class codes")

20,478 polygons → 575,483 of 1,302,622 cells burned with 20 class codes

The grid first, a projected CRS, the burn value

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_30m.tif") as tmpl:
    transform, width, height, crs = tmpl.transform, tmpl.width, tmpl.height, tmpl.crs
print(crs, width, height, transform.a)

EPSG:27700 1214 1073 30.0

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 = 30
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 instead puts the grid's top edge at the data's bottom edge, so the whole grid hangs south of the data and nothing burns.

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

from shapely import box
too_big = rasterize([(box(0, 0, 10, 10), 1000)], out_shape=(10, 10),
                    transform=rasterio.transform.from_origin(0, 10, 1, 1), dtype="uint8")
print("1000 burned into uint8 →", np.unique(too_big))

1 20 int16
1000 burned into uint8 → [255]

A value that does not fit is clamped to the largest the dtype can hold — silently, with no warning. Every code above 255 would come out as 255, merging classes that should be distinct.

4. Understand what happens where shapes overlap

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:

OpenStreetMap land use overlaps a great deal — a nature reserve drawn over a wood, a park over grass — so the order question is real here:

with rasterio.open("elevation_30m.tif") as tmpl:
    grid = {"out_shape": (tmpl.height, tmpl.width), "transform": tmpl.transform, "fill": 0, "dtype": "uint8"}
forward = rasterize(zip(gdf.geometry, gdf["class_code"]), **grid)
backward = rasterize(zip(gdf.geometry[::-1], gdf["class_code"][::-1]), **grid)
print(f"cells whose class depends on row order: {(forward != backward).sum():,} "
      f"of {(forward > 0).sum():,} burned ({(forward != backward).sum() / (forward > 0).sum():.1%})")

cells whose class depends on row order: 23,296 of 575,483 burned (4.0%)

# make priority explicit: burn low priority first, high priority last
priority = {"water": 3, "residential": 2, "farmland": 1, "grass": 0}
gdf["prio"] = gdf["fclass"].map(priority).fillna(-1)
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)
for touched in (True, False):
    road_mask = rasterize(((g, 1) for g in roads.geometry), out_shape=(height, width), transform=transform,
                          fill=0, all_touched=touched, dtype="uint8")
    print(f"all_touched={touched!s:5}  {int(road_mask.sum()):>9,} road cells")

all_touched=True     363,421 road cells
all_touched=False    306,457 road cells

The False figure is not "fewer cells for a thinner line" — the missing cells are gaps, wherever a road passes between cell centres, and a network with gaps is useless for anything that follows the roads. In a dense city most roads cross many cell centres, so the loss is a sixth; for a sparse rural network, or a coarser grid, it is far worse.

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

A study-area mask for Leeds

The data used on this page

Everything is for the Leeds district, in British National Grid:

• elevation_30m.tif — the template grid: the Copernicus 30 m DEM over Leeds, reprojected to 30 m cells.
• leeds.gpkg — the district boundary (ONS, December 2025).
• landcover.gpkg — OpenStreetMap land-use polygons from Geofabrik's 1 January 2026 West Yorkshire extract, with a numeric class_code for each of its fclass values; roads.gpkg and railways.gpkg from the same extract.
• incidents.gpkg — every street-level crime police.uk recorded in Leeds in June 2026.
• layers/ — single-class layers for example 3: woodland, parks, roads, railways, and vineyards (Leeds has a few).

import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.warp import Resampling, calculate_default_transform, reproject
from rasterio.windows import from_bounds

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
leeds = gpd.read_file(f"{ONS}/Local_Authority_Districts_DEC_2025_Boundaries_UK_BFC/FeatureServer/0/query"
                      "?where=LAD25NM%3D%27Leeds%27&outFields=LAD25NM&outSR=27700&f=geojson")
leeds = leeds.set_crs(27700, allow_override=True)
leeds.to_file("leeds.gpkg")
west, south, east, north = leeds.to_crs(4326).total_bounds

DEM = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_W002_00_DEM/"
       "Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif")
with rasterio.open(DEM) as src:
    window = from_bounds(west - 0.02, south - 0.02, east + 0.02, north + 0.02, transform=src.transform)
    window = window.round_offsets().round_lengths()
    data, win_t = src.read(1, window=window), src.window_transform(window)
    t, w, h = calculate_default_transform(src.crs, "EPSG:27700", data.shape[1], data.shape[0],
                                          *rasterio.windows.bounds(window, src.transform), resolution=30)
    dem = np.full((h, w), -9999, dtype="float32")
    reproject(data, dem, src_transform=win_t, src_crs=src.crs, dst_transform=t, dst_crs="EPSG:27700",
              dst_nodata=-9999, resampling=Resampling.bilinear)
with rasterio.open("elevation_30m.tif", "w", driver="GTiff", width=w, height=h, count=1, dtype="float32",
                   crs="EPSG:27700", transform=t, nodata=-9999, compress="deflate") as dst:
    dst.write(dem, 1)

OSM = "https://download.geofabrik.de/europe/united-kingdom/england/west-yorkshire-260101-free.shp.zip"
Path("west-yorkshire.zip").write_bytes(requests.get(OSM, timeout=600).content)
box4326 = (west, south, east, north)
landuse = gpd.read_file("zip://west-yorkshire.zip!gis_osm_landuse_a_free_1.shp", bbox=box4326).to_crs(27700)
landuse = landuse[landuse.intersects(leeds.union_all())].reset_index(drop=True)
landuse["class_code"] = (pd.factorize(landuse["fclass"], sort=True)[0] + 1).astype("uint8")
landuse.to_file("landcover.gpkg")
roads = gpd.read_file("zip://west-yorkshire.zip!gis_osm_roads_free_1.shp", bbox=box4326).to_crs(27700)
roads.to_file("roads.gpkg")
rail = gpd.read_file("zip://west-yorkshire.zip!gis_osm_railways_free_1.shp", bbox=box4326).to_crs(27700)

found = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        found += [(c["id"], float(c["location"]["longitude"]), float(c["location"]["latitude"])) for c in reply.json()]
crimes = pd.DataFrame(found, columns=["id", "lon", "lat"]).drop_duplicates("id")
crimes = gpd.GeoDataFrame(crimes, crs=4326, geometry=gpd.points_from_xy(crimes["lon"], crimes["lat"])).to_crs(27700)
crimes[crimes.within(leeds.union_all())].to_file("incidents.gpkg")

Path("layers").mkdir(exist_ok=True)
landuse[landuse["fclass"] == "forest"].to_file("layers/woodland.gpkg")
landuse[landuse["fclass"] == "park"].to_file("layers/parks.gpkg")
landuse[landuse["fclass"] == "vineyard"].to_file("layers/vineyards.gpkg")
roads.to_file("layers/roads.gpkg")
rail.to_file("layers/railways.gpkg")
print(f"template {w} × {h}; {len(landuse):,} land-use polygons in {landuse['fclass'].nunique()} classes; "
      f"{len(roads):,} roads; {len(gpd.read_file('incidents.gpkg')):,} crimes")

template 1214 × 1073; 20,478 land-use polygons in 20 classes; 116,957 roads; 9,224 crimes

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, blockxsize=256, blockysize=256)     # not the template's strips
    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("leeds.gpkg", "elevation_30m.tif", "leeds_mask.tif")

612,974 of 1,302,622 cells inside (47.1%)

The explicit block size is there because the first version of this function set tiled=True alone and failed with RasterBlockError: The height and width of TIFF dataset blocks must be multiples of 16: the template is a stripped GeoTIFF with blocks one row high, the copied profile kept that block size, and a tile cannot be one row high. 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:

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(f"naive burn: {int(naive.sum()):,}   accumulated: {int(counts.sum()):,}   points: {len(points):,}   "
      f"busiest cell: {int(counts.max())}")

naive burn: 3,968   accumulated: 9,224   points: 9,224   busiest cell: 69

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, blockxsize=256, blockysize=256)
    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

names = build_stack("layers/", "elevation_30m.tif", "predictors.tif")

  ✓ parks.gpkg                  25,536 cells
  ✓ railways.gpkg                8,429 cells  (all_touched)
  ✓ roads.gpkg                 363,421 cells  (all_touched)
  ✓ vineyards.gpkg                  21 cells
  ✓ woodland.gpkg               64,549 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.

Which cells a shape occupies

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.

Fill collisions, dtype overflow and vanishing lines

• 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. In rasterio 1.5 a value too big for the dtype is clamped — 1000 into uint8 gives 255 (step 3) — so distinct large codes merge. Check the maximum against the dtype before burning.

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

Internal links

• The raster data model explained — the transform you must supply
• How to convert a raster to a vector in Python — the reverse conversion, and what it loses
• How to calculate zonal statistics in Python — the same pairing, in the other direction
• How to clip a raster to a polygon in Python — masking built on the same inclusion rules
• How to create a fishnet grid in Python with GeoPandas — the vector equivalent of a grid
• How to choose the right projected CRS for your study area — before rasterizing anything
• How to remove null and empty geometries in GeoPandas — what makes rasterize raise
• How to fix gaps and overlaps in a polygon coverage — when burn order is deciding your answers
• How to measure distance to the grid — rasterising lines for a distance transform

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.