python - Polygonize raster file according to band values?

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
gdal.UseExceptions()

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

    sys.exit(1)

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


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

Answer

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)
        shp.write({
            'geometry': mapping(multipolygon),
            'properties': {'pixelvalue': int(pixel_value)}
        })