Polygonize a Raster to Vector in PyQGIS

Classified rasters — land cover, flood extent, a reclassified slope map, a suitability score in bands — are often needed as polygons: to compute areas per class within districts, to overlay with parcels, to publish as vector data, or to edit by hand. Polygonizing turns each connected group of same-valued pixels into a polygon carrying that value. Done naively on raw output it produces hundreds of thousands of tiny polygons with staircase edges; done with a little preparation it produces a clean, usable layer.

This recipe belongs to Raster Analysis Workflows. It prepares the raster, runs gdal:polygonize, removes NoData and speckle, dissolves by class, smooths edges and computes areas.

Prepare, polygonize, cleanA classified raster with scattered single pixels is first smoothed with a sieve filter that merges regions smaller than a threshold into their neighbours. gdal:polygonize then turns each connected region into a polygon with its class value. The result is cleaned by dropping NoData polygons, dissolving by class if needed, and simplifying or smoothing staircase edges before areas are computed.Fewer, cleaner polygons start in the rasterclassifiedrasterwith specklegdal:sievemerge regions< N pixelsgdal:polygonizeone polygon perconnected regioncleanNoDatadissolvesimplify

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series, with the GDAL provider.
  • A classified integer raster: land-cover codes, flood/no-flood, suitability classes. Continuous rasters such as elevation must be reclassified first, as in reclassifying raster values.

Remove speckle before polygonizing

Classification output is noisy: isolated pixels and tiny clusters that become thousands of meaningless polygons. A sieve filter merges regions smaller than a pixel threshold into their largest neighbour, in the raster, where it is fast.

import processing
from qgis.core import QgsRasterLayer

src = "/data/landcover/classified_2026.tif"
layer = QgsRasterLayer(src, "classified")
px_area = layer.rasterUnitsPerPixelX() * layer.rasterUnitsPerPixelY()
min_area_m2 = 500
threshold = max(1, int(min_area_m2 / px_area))
print(f"pixel {px_area:.0f} m² → sieve threshold {threshold} pixels")

sieved = processing.run("gdal:sieve", {
    "INPUT": src, "THRESHOLD": threshold, "EIGHT_CONNECTEDNESS": False,
    "NO_MASK": False, "MASK_LAYER": None,
    "OUTPUT": "/data/landcover/classified_2026_sieved.tif"})["OUTPUT"]

Breakdown: Expressing the threshold as a minimum area in square metres and converting it to pixels makes the setting meaningful and independent of resolution: 500 m² is five pixels at 10 m resolution and 125 at 2 m. The sieve replaces small regions with the value of their largest neighbour, so it changes the raster — keep the original. Without a sieve, a 10,000 × 10,000 classified image can produce millions of polygons and take hours to polygonize; with one, often a few tens of thousands.

Polygonize

gdal:polygonize walks the raster and builds a polygon for each connected region of equal value, storing the value in a field.

polys = processing.run("gdal:polygonize", {
    "INPUT": sieved, "BAND": 1, "FIELD": "class",
    "EIGHT_CONNECTEDNESS": False, "EXTRA": "",
    "OUTPUT": "/data/landcover/landcover_2026.gpkg"})["OUTPUT"]

from qgis.core import QgsVectorLayer
lc = QgsVectorLayer(polys, "land cover polygons", "ogr")
print(lc.featureCount(), "polygons;", lc.fields().names())

Breakdown: The FIELD name holds the pixel value of each polygon. Writing to GeoPackage rather than shapefile avoids field-name limits and gives a spatial index. Polygon edges follow pixel boundaries exactly, so the output has right-angled "staircase" edges at the raster's resolution — faithful to the data and often ugly; the smoothing section below addresses it. The output CRS is the raster's CRS.

Choose 4- or 8-connectedness

Connectedness decides whether pixels that touch only at a corner belong to the same region. The choice changes both polygon count and shape.

Four versus eight neighboursTwo diagonal pixels of the same class touching only at a corner. With 4-connectedness they are separate regions and become two polygons. With 8-connectedness they are one region, and the resulting polygon touches itself at the shared corner, which some tools treat as invalid. Four-connectedness is the safer default for clean geometry.Do corner-touching pixels connect?4: two polygons8: one, self-touchingdefault: 4valid geometrymore polygons

Four-connectedness (the default) joins pixels only through shared edges; diagonal neighbours stay separate. That produces a few more polygons but always valid geometry. Eight-connectedness joins diagonal neighbours too, which suits thin diagonal features such as streams classified at coarse resolution — but the resulting polygons can touch themselves at a single corner, which strict validators flag as invalid. Use the same setting for the sieve and the polygonize steps so they agree on what a region is.

Drop NoData and unwanted classes

NoData areas become polygons too, unless the raster's NoData value is set — and some classes, such as "unclassified" or "water" in a land analysis, may not be wanted.

from qgis.core import QgsFeatureRequest, edit

drop = lc.getFeatures(QgsFeatureRequest().setFilterExpression('"class" IN (0, 255)'))
ids = [f.id() for f in drop]
with edit(lc):
    lc.deleteFeatures(ids)
print(len(ids), "NoData/unclassified polygons removed;", lc.featureCount(), "remain")

Breakdown: When the raster declares a NoData value, gdal:polygonize skips those pixels automatically; when it does not — common with classification outputs that use 0 or 255 for "nothing" — they become polygons and must be removed. Deleting by class in one call is fast. Doing it before dissolving and simplifying saves work on polygons that would be discarded anyway.

Dissolve and simplify

Polygonize produces one polygon per connected region, so a class appears as many separate polygons. For area statistics per class, dissolving merges them; for display and further vector work, simplifying removes the staircase.

Smoothing staircase edgesA polygon edge following pixel boundaries forms a staircase. Simplification with a tolerance of about half a pixel removes most steps while keeping the shape; a larger tolerance distorts it. Smoothing adds curvature for a natural look but moves edges slightly. Areas should be computed before smoothing if they must match pixel counts.Remove the steps, keep the shapestaircase (pixel edges)simplified (½ pixel)smoothed

dissolved = processing.run("native:dissolve", {
    "INPUT": lc, "FIELD": ["class"], "SEPARATE_DISJOINT": False,
    "OUTPUT": "memory:by_class"})["OUTPUT"]

pixel = layer.rasterUnitsPerPixelX()
simplified = processing.run("native:simplifygeometries", {
    "INPUT": lc, "METHOD": 0, "TOLERANCE": pixel * 0.5,
    "OUTPUT": "memory:landcover_simplified"})["OUTPUT"]

with_area = processing.run("native:fieldcalculator", {
    "INPUT": simplified, "FIELD_NAME": "area_ha", "FIELD_TYPE": 0,
    "FIELD_LENGTH": 12, "FIELD_PRECISION": 3, "FORMULA": "$area / 10000",
    "OUTPUT": "memory:landcover_final"})["OUTPUT"]

Breakdown: Dissolving by class gives one (multi)polygon per class — the right input for "how many hectares of forest", but too coarse for most mapping. Simplifying with a tolerance of half a pixel removes the steps without moving edges by more than the data's own resolution; much larger tolerances visibly distort shapes. native:smoothgeometry gives a softer look for display but moves edges further. Simplifying polygons independently can open small gaps between neighbours; for a strict coverage, simplify with topology-aware tools or accept the staircase. Simplifying geometry compares the methods.

Summarise classes within districts

The usual reason to polygonize is to combine the classes with other vector data: hectares of each land-cover class per municipality, flooded area per parcel. Once the classes are polygons, an overlay and a group-by answer it.

districts = QgsVectorLayer("/data/admin/districts.gpkg|layername=districts", "districts", "ogr")
pieces = processing.run("native:intersection", {
    "INPUT": with_area, "OVERLAY": districts,
    "INPUT_FIELDS": ["class"], "OVERLAY_FIELDS": ["district"],
    "OUTPUT": "memory:"})["OUTPUT"]

table = defaultdict(float)
for f in pieces.getFeatures():
    table[(f["district"], f["class"])] += f.geometry().area() / 1e4
for (district, cls), ha in sorted(table.items())[:12]:
    print(f"{district:<20} class {cls:>3}  {ha:8.1f} ha")

Breakdown: Intersecting the class polygons with district boundaries splits each polygon where it crosses a boundary and keeps one attribute from each side. Recomputing area on the pieces, rather than reusing the area_ha field from before the overlay, gives the area inside each district. For this particular question there is also a raster-only route — zonal statistics with a categorical histogram counts pixels per class per district without polygonizing at all, and is faster on large rasters. Polygonizing pays off when the polygons themselves are needed afterwards: for overlay with parcels, for editing, or for publication.

Check areas against pixel counts

Polygon areas should agree with pixel counts per class. Comparing them catches lost polygons, unexpected NoData and distortion from simplification.

from collections import defaultdict

area_by_class = defaultdict(float)
for f in lc.getFeatures():
    area_by_class[f["class"]] += f.geometry().area()

stats = processing.run("native:rasterlayeruniquevaluesreport", {
    "INPUT": sieved, "BAND": 1, "OUTPUT_TABLE": "memory:"})["OUTPUT_TABLE"]
for row in stats.getFeatures():
    cls, m2 = row["value"], row["count"] * px_area
    poly = area_by_class.get(cls, 0)
    print(f"class {cls:>3}: raster {m2 / 1e4:10.1f} ha  polygons {poly / 1e4:10.1f} ha")

Breakdown: native:rasterlayeruniquevaluesreport counts pixels per value; multiplying by the pixel area computed earlier converts counts to square metres. Before simplification, polygon areas must match exactly; after simplification they differ slightly, and the comparison shows by how much. A class missing from the polygons usually means it was deleted as NoData by mistake.

QGIS version compatibility

gdal:sieve, gdal:polygonize and the native vector algorithms work on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. SEPARATE_DISJOINT on native:dissolve was added in 3.32; omit it on older releases. native:rasterlayeruniquevaluesreport has existed since 3.4.

Troubleshooting

  • Millions of polygons. Sieve first, with a threshold tied to a minimum mapping area.
  • A huge polygon covers the whole extent. NoData was not declared; delete the NoData-class polygon or set NoData on the raster.
  • Invalid geometries after polygonizing. 8-connectedness created self-touching rings; use 4-connectedness or fix geometries.
  • Gaps between polygons after simplifying. Polygons were simplified independently; use a smaller tolerance.

Conclusion

Sieve speckle away in the raster with a threshold from a minimum area, polygonize with 4-connectedness into GeoPackage, drop NoData and unwanted classes, dissolve for per-class statistics, simplify by about half a pixel for clean edges, and check areas against pixel counts.

Frequently Asked Questions

Can I polygonize a continuous raster? Not usefully — every distinct value becomes a polygon. Reclassify into bands first.

How do I polygonize only one class? Reclassify everything else to NoData, then polygonize; or polygonize all and filter.

Is there a faster way for huge rasters? Process in tiles with a small overlap and merge results, or use gdal_polygonize.py directly with a large cache.

Should I keep the raster after polygonizing? Yes. The raster is the source of truth for the classification; polygons are a derived, simplified view. Keep both, and record the sieve threshold and simplification tolerance with the polygons so the derivation can be repeated.

Can I attach the class names to the polygons? Join a lookup table of class codes and names on the class field, so the layer is readable without the classification's documentation.

Why do polygons look blocky at high zoom? They follow pixel edges. Simplify or smooth for display; keep the originals for analysis.