0

I want to polygonize a raster by using GDAL in python. I found this webpage below: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#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 multi polygon)

EDIT Here ıs my polygonıze 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 dıfferent band valuesi I should get a layer wıth 500 features

  • You can union the separate polygons having the same value afterwards. – user30184 Sep 9 '18 at 20:59
  • How can I do that after polygonize? – Mustafa Uçar Sep 9 '18 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. – Michael Stimson Sep 10 '18 at 3:56
  • @MichaelStimson I add my polygonıze code – Mustafa Uçar Sep 10 '18 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 – Michael Stimson Sep 10 '18 at 5:54
2

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.

Code

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)}
        })
  • And how can I add more fields in attribute table? – Mustafa Uçar Sep 21 '18 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 '18 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 – Mustafa Uçar Sep 25 '18 at 4:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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