gdal_polygonize creating empty line file instead of polygon

Question

I'm not sure what I'm doing wrong here. I suspect it is something simple and a duh moment. I'm just trying to get gdal_polygonize.py to make a polygon shapefile of a raster with a single value. It creates a SHP, but it is an empty line file.

Thoughts?

from osgeo import gdal, ogr
import os
import subprocess
"""
# Utility function to polygonize a raster with only 1 value
# @param
# [in_raster] - full path to a GeoTIFF
# Example: rast = r'C:\Temp\MyRast.tif'
# [shp] - full path to a shapefile to be created
# Example: shp = r'C:\Temp\MyShape.shp'
# [region] - integer indicating which connectedness setting to use
# Example call: polygonize(r'C:\Temp\MyRast.tif', r'C:\Temp\MyShape.shp', 4)
"""
def polygonize (in_raster, shp, region):
 ## Open the raster as a GDAL object
 raster = gdal.Open(in_raster)
 ## Get the projection of the raster
 prj = raster.GetProjection()
 ## Create a spatial reference system object from the projection
 shp_srs = osr.SpatialReference(wkt=prj)
 ## Set the driver for writing the SHP
 driver = ogr.GetDriverByName('ESRI Shapefile')
 ## Specify the output SHP dataset
 out_dataset = driver.CreateDataSource(shp)
 ## Create a layer object for GDAL to write to with the PRJ set the same as input raster
 out_layer = out_dataset.CreateLayer(os.path.basename(shp)[:-4], srs=shp_srs)
 if region == 4:
 ## Create a variable to hold the command syntax for 'gdal_polygonize' operation
 cmd = ['python', r'c:\Program Files\QGIS 2.18\bin\gdal_polygonize.py', in_raster, '-f', 'ESRI Shapefile', shp, out_layer, 'DN']
 ## Execute the command using subprocess
 process = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
 ## Wait for the process to finish 
 out, err = process.communicate()
 ## Let user know what is happening
 print ("Finished polygonizing: {0}".format(os.path.basename(shp)))
 ## Get error code
 errcode = process.returncode
 ## Print error information if wanted
 #print (out, err)
 ## Return the error code if needed - non-zero is error
 #return errcode

Here is an example input raster.

• dgj32784: So why would a call to gdal.Polygonize that looks like this result in a shapefile with geometry for all values and the code above only providing geometry where it should? gdal.Polygonize( src_band, None, out_layer, -1, [], None ) Obviously preface this with selecting the band from the raster. src_band = raster.GetRasterBand(1)

Answer 1

Using the QGIS GUI, it works for me with this command line:

gdal_polygonize.bat F:/path/to/MyRast/MyRast.tif -f "ESRI Shapefile" G:/test.shp test DN

Your codes misses the output field name which defaults to DN in the GUI.

• dgj32784: Thanks Andre, your post made me reevaluate what I was passing. Sloppy work on my part. I had copy/pasted a variable assignment and simply added the new name to the left without deleting the old. So the raster I thought I was passing in wasn't actually a raster. That's the main reason it wasn't doing anything. But it does appear that 'DN' is required in my case for the polygonization to occur, which I'm a little confused about since after adding it my output file name is 'DN.shp'. What's causing that?
• AndreJ: You have still mixed up out_layer and dst_layer.
• dgj32784: Right, that was a copy/paste issue on the post. I changed variable names. So I'm getting an empty shapefile made that has the name I want and a shapefile named 'DN.shp' that actually has the polygons. Any idea what's going on? I'm not sure I have the command right...
• AndreJ: I don't have much experience with python, but do you tell the driver to create a polygon layer anywhere?
• dgj32784: Driver and related shp/layer writing should be handled with this portion: ## Set the driver for writing the SHP driver = ogr.GetDriverByName('ESRI Shapefile') ## Specify the output SHP dataset out_dataset = driver.CreateDataSource(shp) ## Create a layer object for GDAL to write to with the PRJ set the same as input raster out_layer = out_dataset.CreateLayer(os.path.basename(shp)[:-4], srs=shp_srs)