17

I am having trouble using the Polygonize function in python. The cookbook example for this can be found here.

The relevant portion of my code is:

sourceRaster = gdal.Open('myraster.tif')
band = sourceRaster.GetRasterBand(1)
bandArray = band.ReadAsArray()
outShapefile = "polygonized"
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outShapefile+".shp"):
    driver.DeleteDataSource(outShapefile+".shp")
outDatasource = driver.CreateDataSource(outShapefile+ ".shp")
outLayer = outDatasource.CreateLayer("polygonized", srs=None)
gdal.Polygonize( band, None, outLayer, -1, [], callback=None )
outDatasource.Destroy()
sourceRaster = None

I know that the band has relevant information, here is a snippet of bandArray:

array([[ 4.,  4.,  3.,  3.,  3.,  2.,  2.,  2.,  2.,  3.,  3.,  3.,  3.,
         3.,  3.,  3.,  3.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,
         4.,  4.,  4.,  4.],

When I open the attribute table in QGIS, it is empty: QGIS screencapture

Edit:

The conversion works just fine in QGIS using Raster -> Conversion -> Polygonize tool

Screenshot of raster to be polygonized:

raster to polygonize

And screenshot of resultant conversion from QGIS tool:

polygonized raster from QGIS tool

I am using the Enthought distro on Windows 7, GDAL version 1.10.0-3

The problem is that I cannot polygonize a raster in python using GDAL and the cookbook example, I can polygonize this same raster with no problem in the QGIS GUI

6
  • What does your raster look like? Does it really contain polygons? Does it work if you use gdal_polygonize.py instead?
    – BradHards
    Apr 7, 2014 at 11:24
  • Edited to add screenshots of working process in QGIS Apr 7, 2014 at 11:57
  • What's the actual problem here?
    – Fezter
    Apr 7, 2014 at 12:06
  • Added specific problem Apr 7, 2014 at 12:14
  • 4
    I had a similar problem (blank shapefile being created), and creating the field didn't help. What I was doing wrong was that I hadn't closed the shapefile in my code before calling polygonize. You do close it in your example, I am just posting this for others' reference.
    – Stephanie
    Nov 6, 2015 at 21:16

2 Answers 2

23

The problem is that I was not creating a field to store the raster band. After digging through the gdal_polygonize.py file, I realized this is not automatically done when calling gdal.Polygonize, which instead uses the function found here.

Here is the extra step needed to create a field and write a band to the field:

newField = ogr.FieldDefn('MYFLD', ogr.OFTInteger)
outLayer.CreateField(newField)

We can then write the band to this field, with an index of 0:

gdal.Polygonize(band, None, outLayer, 0, [], callback=None )
4
  • I am also trying to use gdal.Polygonize() function to get my raster as polygon in python. But at the end line it is showing runtime error!! why? Jul 20, 2016 at 9:04
  • It works well with georeferenced raster files. The result is too many polygons, but i only want one big polygon showing the outline of raster. Does anyone have any idea how it works at the same time? Jul 20, 2016 at 9:58
  • Still i'm getting empty shapefile, but i have rows in dbf file. Please clearify me! Jun 6, 2017 at 10:46
  • I just had this problem but instead of adding a dummy field, you can enter an index of -1. See here that field is only added if index >=0.
    – jon_two
    Jun 22, 2017 at 20:40
4

I had the same problem and I fixed by closing the shapefile after using poligonize ("dst_ds.Destroy())":

def polygon_response(raster, poligonized_shp):
src_ds = gdal.Open(raster)
if src_ds is None:
    print('Unable to open %s' % src_filename)
    sys.exit(1)

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

#  create output datasource
dst_layername = poligonized_shp
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 )
dst_ds.Destroy()

polygon_response(raster, poligonized_shp)
2
  • this worked for me, thanks a lot! Feb 10, 2021 at 15:44
  • this just saved me. Thanks Oct 20, 2021 at 20:46

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.