I want to polygonize a raster by using GDAL in python. I found this webpage: Polygonize a Raster Band. But it creates a polygon feature layer. I want every band with the same value to create one feature.

For instance, even if band values are far away from each other but the same value, they should be one feature (maybe multipolygon)

EDIT Here is my polygonize code:

from osgeo import gdal, ogr
import sys
# this allows GDAL to throw Python Exceptions

#  get raster datasource
src_ds = gdal.Open( "INPUT.tif" )
if src_ds is None:
    print 'Unable to open %s' % src_filename

    srcband = src_ds.GetRasterBand(3)
except RuntimeError, e:
    # for example, try GetRasterBand(10)
    print 'Band ( %i ) not found' % band_num
    print e

#  create output datasource
dst_layername = "POLYGONIZED_STUFF"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

If there are 500 different band values I should get a layer with 500 features.

  • You can union the separate polygons having the same value afterwards.
    – user30184
    Sep 9, 2018 at 20:59
  • How can I do that after polygonize? Sep 9, 2018 at 21:14
  • It depends what software you're using other than GDAL. I can assume you're using GDAL in python but without any code to base an opinion on I can't make a suggestion on how to dissolve all the parts into a multipolygon; I can't even be certain your geometry storage supports multipolygon features. Can you post a snippet of your code please. Sep 10, 2018 at 3:56
  • @MichaelStimson I add my polygonıze code Sep 10, 2018 at 5:45
  • This might help gis.stackexchange.com/questions/149959/… or for something a little simpler darrencope.com/2015/01/23/dissolve-shapefiles-using-ogr Sep 10, 2018 at 5:54

1 Answer 1


It could be easily done with Rasterio, Shapely, and Fiona ; using the following processing scheme:

  1. Read input band with Rasterio
  2. Get unique values of the input band
  3. Polygonize input band with Rasterio
  4. Get a list of polygons for each unique value...
  5. ...which can be converted into a MultiPolygon with shapely
  6. Write each record to an output shapefile with Fiona.


import numpy as np
import fiona
import rasterio
import rasterio.features
from shapely.geometry import shape, mapping
from shapely.geometry.multipolygon import MultiPolygon

# Read input band with Rasterio
with rasterio.open('INPUT.tif') as src:
    crs = src.crs
    src_band = src.read(3)
    # Keep track of unique pixel values in the input band
    unique_values = np.unique(src_band)
    # Polygonize with Rasterio. `shapes()` returns an iterable
    # of (geom, value) as tuples
    shapes = list(rasterio.features.shapes(src_band, transform=src.transform))

shp_schema = {
    'geometry': 'MultiPolygon',
    'properties': {'pixelvalue': 'int'}

# Get a list of all polygons for a given pixel value
# and create a MultiPolygon geometry with shapely.
# Then write the record to an output shapefile with fiona.
# We make use of the `shape()` and `mapping()` functions from
# shapely to translate between the GeoJSON-like dict format
# and the shapely geometry type.
with fiona.open('output.shp', 'w', 'ESRI Shapefile', shp_schema, crs) as shp:
    for pixel_value in unique_values:
        polygons = [shape(geom) for geom, value in shapes
                    if value == pixel_value]
        multipolygon = MultiPolygon(polygons)
            'geometry': mapping(multipolygon),
            'properties': {'pixelvalue': int(pixel_value)}
  • And how can I add more fields in attribute table? Sep 21, 2018 at 9:31
  • If you want to add more field to your output shapefile, you have to: 1. Update the shapefile schema shp_schema['properties'] with your fields ; e.g. shp_schema['properties']['my_field'] = 'str' 2. Assign the value in the for loop, e.g. like I dit it for pixelvalue.
    – Y. Forget
    Sep 23, 2018 at 13:04
  • Actually I don't know how many attribute user add. So I can only keep them in list and I should add these attributes by iterating. Can you give me an example? I also have a database table which holds pixel value and other attributes values corresponding to pixel values, I must also add them by using this Sep 25, 2018 at 4:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.