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

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
camdenl
  • 1,213
  • 2
  • 12
  • 27

2 Answers2

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 )
camdenl
  • 1,213
  • 2
  • 12
  • 27
  • 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? – Shiuli Pervin Jul 20 '16 at 09: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? – Shiuli Pervin Jul 20 '16 at 09:58
  • Still i'm getting empty shapefile, but i have rows in dbf file. Please clearify me! – Satya Chandra Jun 06 '17 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 '17 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)
Fjord
  • 41
  • 4