How to Convert a Raster to a Vector in Python

On this page

Problem statement

Polygonising a raster turns every run of same-valued pixels into a polygon — which is exactly right for a classified land cover grid or a water mask, and produces a staircase of tiny shapes if you skip the cleaning. It is the step behind:

extracting land cover classes from a classified TIFF
converting a water or flood mask into polygons
turning valid raster regions into a GeoPackage, Shapefile or GeoJSON
creating vector boundaries from a suitability or segmentation raster

In Python, the standard workflow uses:

rasterio to read the raster
rasterio.features.shapes to polygonize pixel regions
GeoPandas to store and export the result
Shapely to build geometry objects

This works best for categorical rasters where neighbouring pixels share the same class value. Every example below runs on a real one: ESA WorldCover's 10 m land cover for 2021, cut to Rutland Water and Oakham.

Quick answer

Polygonizing a raster merges neighbouring cells that share a value into vector polygons.

To convert a raster to vector polygons in Python, the usual workflow is:

open the raster with Rasterio
read the raster band and define a valid-data mask
extract polygons with rasterio.features.shapes
load the results into a GeoDataFrame
save to GeoPackage, Shapefile or GeoJSON

Basic example:

import os
import rasterio
from rasterio.features import shapes
import geopandas as gpd
from shapely.geometry import shape

raster_path = "data/worldcover_rutland.tif"
output_path = "output/landcover_polygons.gpkg"

os.makedirs("output", exist_ok=True)

with rasterio.open(raster_path) as src:
    band = src.read(1)
    crs = src.crs
    mask = band != src.nodata if src.nodata is not None else band > 0

    records = []
    for geom, value in shapes(band, mask=mask, transform=src.transform):
        records.append({"geometry": shape(geom), "class_value": int(value)})

gdf = gpd.GeoDataFrame(records, crs=crs)
gdf.to_file(output_path)
print(len(gdf), "polygons from", band.size, "pixels; CRS", gdf.crs.to_string())

5003 polygons from 1958400 pixels; CRS EPSG:4326

This approach is best for classified rasters. For continuous rasters like elevation or imagery, reclassify first or you may create too many polygons.

Inspect values, drop background, polygonise, keep class meanings

Skipping the mask step is what produces one enormous nodata polygon.

Step 1: Load the raster and inspect its values

Before polygonizing, check the raster metadata and values. You need to know:

CRS
affine transform
nodata value
raster class values

import numpy as np

with rasterio.open(raster_path) as src:
    band = src.read(1)
    print("CRS:", src.crs)
    print("Transform:", tuple(round(v, 6) for v in src.transform[:6]))
    print("NoData:", src.nodata, "| dtype:", src.dtypes[0], "| shape:", band.shape)

    values, counts = np.unique(band, return_counts=True)
    print({int(v): int(c) for v, c in zip(values, counts)})

CRS: EPSG:4326
Transform: (8.3e-05, 0.0, -0.72, 0.0, -8.3e-05, 52.7)
NoData: 0.0 | dtype: uint8 | shape: (960, 2040)
{10: 317017, 30: 518363, 40: 835595, 50: 45976, 60: 12480, 80: 225992, 90: 2977}

The raster has only seven values — WorldCover's class codes, where 10 is tree cover, 30 grassland, 40 cropland, 50 built-up, 60 bare, 80 permanent water and 90 wetland — so it is suitable for polygonizing. If it had many unique values, it would probably be continuous data and should be reclassified first.

For very large rasters, np.unique() can be slow or memory-intensive. In that case, inspect a smaller window first or read the class definitions from the product's documentation.

Step 2: Filter out nodata or background pixels

You usually do not want nodata or background cells turned into polygons. Build a mask so only valid cells are included.

with rasterio.open(raster_path) as src:
    band = src.read(1)
    if src.nodata is not None:
        mask = band != src.nodata
    else:
        mask = band > 0  # example if 0 is background

print(f"{int(mask.sum()):,} of {mask.size:,} pixels are valid")

1,958,400 of 1,958,400 pixels are valid

The mask controls which cells are polygonized. This cut contains no nodata at all; at the edge of a tile, or around the sea, the mask is what stops the empty cells becoming one enormous polygon.

Step 3: Convert raster regions to vector geometries

Use rasterio.features.shapes to extract polygons. It returns geometry/value pairs. Adjacent pixels with the same value are grouped into polygon features.

with rasterio.open(raster_path) as src:
    band = src.read(1)
    mask = band != src.nodata if src.nodata is not None else band > 0

    results = shapes(band, mask=mask, transform=src.transform)
    geom, value = next(iter(results))
    print(value, geom["type"], geom["coordinates"][0][:2])

30.0 Polygon [(-0.6559999999999997, 52.7), (-0.6559999999999997, 52.69991666666667)]

The transform is required so raster row and column positions become real map coordinates — here longitude and latitude, because WorldCover is delivered in EPSG:4326.

Step 4: Build a GeoDataFrame from the extracted shapes

Convert the GeoJSON-like geometry dictionaries into Shapely geometries, then create a GeoDataFrame.

with rasterio.open(raster_path) as src:
    band = src.read(1)
    crs = src.crs
    mask = band != src.nodata if src.nodata is not None else band > 0

    records = []
    for geom, value in shapes(band, mask=mask, transform=src.transform):
        records.append({
            "geometry": shape(geom),
            "class_value": int(value)
        })

gdf = gpd.GeoDataFrame(records, crs=crs)
print(gdf.head(3))
print(gdf["class_value"].value_counts().sort_index().to_dict())

                                            geometry  class_value
0  POLYGON ((-0.656 52.7, -0.656 52.69992, -0.655...           30
1  POLYGON ((-0.65292 52.7, -0.65292 52.69992, -0...           30
2  POLYGON ((-0.63217 52.7, -0.63217 52.69992, -0...           40
{10: 1169, 30: 2067, 40: 440, 50: 926, 60: 208, 80: 104, 90: 89}

Now each polygon has a class_value attribute from the raster.

Step 5: Save the output as GeoPackage, Shapefile or GeoJSON

Export the GeoDataFrame in the format you need.

gdf.to_file("output/landcover_polygons.gpkg", layer="landcover")
gdf.to_file("output/landcover_polygons.shp")
gdf.to_file("output/landcover_polygons.geojson", driver="GeoJSON")
for path in ("output/landcover_polygons.gpkg", "output/landcover_polygons.shp", "output/landcover_polygons.geojson"):
    print(f"{path:36} {os.path.getsize(path) / 1e6:5.1f} MB")

output/landcover_polygons.gpkg         8.0 MB
output/landcover_polygons.shp          3.3 MB
output/landcover_polygons.geojson      7.4 MB

GeoPackage is one file and keeps field names intact; a shapefile is widely supported but is really several files — the .shp measured above is only the geometry, with the attributes in a .dbf beside it; GeoJSON is often easier for web mapping and data exchange.

Polygonising WorldCover classes

Examples 1 and 3 in one picture: polygonizing is the big reduction, and the area filter removes most of what is left.

The data used on this page

ESA WorldCover 10 m 2021 v200 (© ESA WorldCover project 2021 / Contains modified Copernicus Sentinel data (2021) processed by ESA WorldCover consortium, CC BY 4.0), read straight from its cloud-optimised GeoTIFF on AWS for a 0.17° × 0.08° window over Rutland Water and Oakham — 2,040 × 960 pixels — and, for the continuous-raster comparison, the Copernicus GLO-30 DEM over the same window.

import os
import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize, shapes
from rasterio.windows import from_bounds
from shapely.geometry import shape

Path("data").mkdir(exist_ok=True)
BOUNDS = (-0.72, 52.62, -0.55, 52.70)          # Rutland Water and Oakham, lon/lat
SOURCES = {
    "data/worldcover_rutland.tif": "https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
                                   "ESA_WorldCover_10m_2021_v200_N51W003_Map.tif",
    "data/dem_rutland.tif": "https://copernicus-dem-30m.s3.amazonaws.com/"
                            "Copernicus_DSM_COG_10_N52_00_W001_00_DEM/Copernicus_DSM_COG_10_N52_00_W001_00_DEM.tif",
}
for path, url in SOURCES.items():
    if not Path(path).exists():
        with rasterio.open(url) as src:
            window = from_bounds(*BOUNDS, src.transform)
            profile = src.profile | {"driver": "GTiff", "height": int(round(window.height)),
                                     "width": int(round(window.width)), "transform": src.window_transform(window),
                                     "tiled": False, "compress": "deflate"}
            profile.pop("blockxsize", None), profile.pop("blockysize", None)
            with rasterio.open(path, "w", **profile) as dst:
                dst.write(src.read(1, window=window), 1)
print("rasterio", rasterio.__version__, "| GDAL", rasterio.__gdal_version__, "|",
      {p: rasterio.open(p).shape for p in SOURCES})

rasterio 1.5.1 | GDAL 3.12.4 | {'data/worldcover_rutland.tif': (960, 2040), 'data/dem_rutland.tif': (288, 408)}

Example 1: Convert a classified raster, and keep what the classes mean

CLASSES = {10: "tree cover", 20: "shrubland", 30: "grassland", 40: "cropland", 50: "built-up",
           60: "bare", 70: "snow and ice", 80: "permanent water", 90: "wetland", 95: "mangroves", 100: "moss"}

with rasterio.open("data/worldcover_rutland.tif") as src:
    band = src.read(1)
    mask = band != src.nodata if src.nodata is not None else band > 0
    features = [
        {"geometry": shape(geom), "class_value": int(value), "class_name": CLASSES[int(value)]}
        for geom, value in shapes(band, mask=mask, transform=src.transform)
    ]
    crs = src.crs

landcover = gpd.GeoDataFrame(features, crs=crs).to_crs(27700)
landcover["area_ha"] = landcover.area / 10_000
summary = landcover.groupby("class_name").agg(polygons=("area_ha", "size"), hectares=("area_ha", "sum"))
print(summary.sort_values("hectares", ascending=False).round(1).to_string())

                 polygons  hectares
class_name                         
cropland              440    4366.6
grassland            2067    2709.2
tree cover           1169    1656.6
permanent water       104    1181.2
built-up              926     240.3
bare                  208      65.2
wetland                89      15.6

Reprojecting to British National Grid before measuring is what makes the hectares real; in the raster's own CRS, .area would be in square degrees.

Example 2: Convert only one raster class to vector

If you only want one value, such as the water class, create a boolean mask for that class.

target_value = 80                                 # permanent water

with rasterio.open("data/worldcover_rutland.tif") as src:
    band = src.read(1)
    class_mask = band == target_value

    water = gpd.GeoDataFrame(
        [{"geometry": shape(geom), "class_value": int(value)}
         for geom, value in shapes(band, mask=class_mask, transform=src.transform)],
        crs=src.crs,
    ).to_crs(27700)

water["area_ha"] = water.area / 10_000
largest = water.sort_values("area_ha", ascending=False).iloc[0]
print(f"{len(water)} water polygons; the largest is {largest.area_ha:,.0f} ha with "
      f"{len(largest.geometry.interiors)} holes; {int((water.area_ha < 0.1).sum())} are under 0.1 ha")
water.to_file("output/water.gpkg", layer="water")

104 water polygons; the largest is 1,142 ha with 57 holes; 78 are under 0.1 ha

With a mask, the value check inside the loop that you often see is redundant: shapes only yields regions under the mask, and here that is class 80. The largest polygon is Rutland Water, and its holes are islands and patches of shoreline the classifier called something other than water.

Example 3: Dissolve by class and drop specks

Polygonized output is noisy: every isolated pixel is its own polygon.

area_threshold_ha = 0.5
kept = landcover[landcover["area_ha"] >= area_threshold_ha]
dissolved = kept.dissolve(by="class_name", aggfunc={"area_ha": "sum"})
print(f"{len(landcover):,} polygons -> {len(kept):,} of at least {area_threshold_ha} ha "
      f"({len(landcover) - len(kept):,} specks, {landcover.loc[landcover.area_ha < area_threshold_ha, 'area_ha'].sum():.0f} ha) "
      f"-> {len(dissolved)} dissolved features")
print("invalid geometries after polygonizing:", int((~landcover.is_valid).sum()))

5,003 polygons -> 601 of at least 0.5 ha (4,402 specks, 244 ha) -> 7 dissolved features
invalid geometries after polygonizing: 0

Filtering by area must happen in a metric CRS, which is why landcover was reprojected in example 1.

Example 4: 4-connectivity, 8-connectivity, and dtypes

with rasterio.open("data/worldcover_rutland.tif") as src:
    band, transform = src.read(1), src.transform

for connectivity in (4, 8):
    n = sum(1 for _ in shapes(band, connectivity=connectivity, transform=transform))
    print(f"connectivity={connectivity}: {n:,} polygons")

for dtype in ("uint8", "int64", "float64", "bool"):
    try:
        n = sum(1 for _ in shapes(band.astype(dtype), transform=transform))
        print(f"{dtype:8} accepted, {n:,} polygons")
    except ValueError as exc:
        print(f"{dtype:8} ValueError: {str(exc)[:70]}")

connectivity=4: 5,003 polygons
connectivity=8: 4,146 polygons
uint8    accepted, 5,003 polygons
int64    accepted, 5,003 polygons
float64  accepted, 5,003 polygons
bool     ValueError: image dtype must be one of: int8, int16, int32, int64, uint8, uint16,

shapes merges only pixels that share an edge by default; connectivity=8 also merges pixels that touch at a corner, which here joined about a sixth of the polygons into their diagonal neighbours. On dtypes, this rasterio accepted 64-bit integers and floats — older releases did not — and refused a boolean array outright; cast a mask to uint8 before polygonizing it.

Example 5: A continuous raster, and the round trip back

with rasterio.open("data/dem_rutland.tif") as dem:
    elevation = dem.read(1)
    start = time.perf_counter()
    raw = sum(1 for _ in shapes(elevation, transform=dem.transform))
    raw_seconds = time.perf_counter() - start
    bands = (np.floor(elevation / 25) * 25).astype("int16")          # 25 m elevation bands
    banded = sum(1 for _ in shapes(bands, transform=dem.transform))
print(f"DEM, {elevation.size:,} cells: {raw:,} polygons raw ({raw_seconds:.1f} s), {banded:,} after 25 m banding")

with rasterio.open("data/worldcover_rutland.tif") as src:
    burned = rasterize(zip(gpd.GeoDataFrame(features, crs=src.crs).geometry,
                           [f["class_value"] for f in features]),
                       out_shape=src.shape, transform=src.transform, dtype="uint8")
    print("rasterize(shapes(band)) == band:", bool((burned == src.read(1)).all()))

DEM, 117,504 cells: 103,698 polygons raw (0.3 s), 296 after 25 m banding
rasterize(shapes(band)) == band: True

Polygonizing a continuous surface gives nearly a polygon per cell, which is why reclassifying first matters. And rasterize is the inverse: burning the polygons back into the same grid reproduced every one of the original pixels.

Adjacent equal pixels become polygons

When you polygonize a raster in Python, the raster is scanned for groups of adjacent pixels with the same value. Each region becomes a polygon.

The pixel value becomes an attribute on the polygon it produced.

The key parts are:

raster band values: define which class each cell belongs to
mask: limits extraction to valid cells
transform: converts pixel coordinates into map coordinates
CRS: keeps the output aligned with other GIS layers

rasterio.features.shapes does not create one polygon per cell unless every cell is isolated. It groups neighbouring cells with the same value into larger polygons — 5,003 of them from nearly two million cells of land cover here.

This is different from other raster-to-vector tasks:

polygonizing creates area features
point conversion creates a point for each cell
contour extraction creates lines from continuous surfaces

For most workflows, polygonizing is the right method when the raster represents zones or classes.

Continuous rasters, noise and nodata

Continuous rasters can produce too many polygons

Elevation rasters, temperature grids, and imagery contain many unique values. Polygonized directly, the DEM in example 5 produced nearly a polygon per cell. Reclassify first.

Small pixel regions create noisy output

Classified rasters contain isolated cells and slivers, and each becomes a small polygon. A common follow-up, as in example 3, is to:

filter small areas in a metric CRS
dissolve polygons by class
smooth boundaries if appropriate

Nodata handling affects the result

If nodata is missing or incorrect, you may get unwanted polygons around empty cells. Always check src.nodata and confirm the mask logic. This WorldCover cut reported a nodata value of 0 and contained none; the band > 0 fallback in the mask line is for rasters that declare no nodata at all.

CRS must be preserved

The output GeoDataFrame should use the raster CRS:

gdf = gpd.GeoDataFrame(records, crs=crs)

If the output does not line up with other layers, the problem is usually CRS assignment or later reprojection.

Invalid geometries can appear

Polygonized output can include invalid geometries on complex raster edges; this land-cover cut produced none, as example 3 checked. Check validity if later overlay operations fail:

invalid = ~gdf.is_valid
print(gdf[invalid])

Large rasters can be slow to polygonize

Very large rasters can be slow and memory-intensive to process, and the output can contain complex geometries with many holes — the reservoir polygon in example 2 has 57 holes. Clip the raster to the area of interest first, as this page's setup does with a window read.

Internal links

If you need background on when this workflow makes sense, see Raster vs vector data in GIS.

Related next steps:

If your output layer does not line up with other data, see How to fix CRS mismatch errors in GeoPandas.

FAQ

How do I convert a raster to polygons in Python?

Use rasterio to read the raster and rasterio.features.shapes to extract polygon geometries from grouped pixel regions. Then store the result in a GeoDataFrame and export it.

What Python library is used to polygonize a raster?

The standard tool is rasterio.features.shapes. GeoPandas is commonly used after that to manage and save the vector output.

Can I convert only one raster value to vector polygons?

Yes. Build a mask such as band == 1 and polygonize only that class. This is common for binary masks like flooded versus not flooded.

Why does raster to vector conversion create too many polygons?

This usually happens when the raster is continuous or noisy. Many unique values or isolated pixels create many separate polygon regions. Reclassification or filtering is often needed first.

How do I save polygonized raster output as a Shapefile or GeoJSON?

Use GeoPandas:

gdf.to_file("output.shp")
gdf.to_file("output.geojson", driver="GeoJSON")

Shapefile is widely supported, while GeoJSON is often easier for web and exchange workflows.

How do I dissolve the polygonized output by class value?

After building the GeoDataFrame, group by the class attribute and merge geometries: gdf.dissolve(by="class_value"). This combines all polygons of the same class into one multipolygon feature.

How do I remove tiny sliver polygons after polygonizing?

Compute area in a projected CRS and filter by a threshold, for example gdf = gdf[gdf.geometry.area > 100]. Make sure the layer is in a metric CRS such as a UTM zone first, or the area values will be in degrees.

Does `rasterio.features.shapes` use 4-connectivity or 8-connectivity?

By default it uses 4-connectivity, so only pixels sharing an edge are merged into the same polygon. Pass connectivity=8 to also merge pixels that touch only at a corner — on the WorldCover window used here that cut 5,003 polygons to 4,146.

Why does `rasterio.features.shapes` raise a dtype error on my band?

shapes accepts only numeric dtypes. rasterio 1.5 took uint8, int64 and float64 alike and refused a boolean array — the usual culprit is a mask from a comparison such as band == 80, which needs .astype("uint8"). Older rasterio releases also refused 64-bit types; cast to int32 or float32 if yours does.

How do I go the other way and rasterize vector polygons back to a raster?

Use rasterio.features.rasterize, passing (geometry, value) pairs with an out_shape and transform. It is the inverse of shapes and burns polygon values into a pixel grid.