2

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.

  • 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) – dgj32784 Mar 12 '17 at 4:38
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.

  • 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? – dgj32784 Mar 11 '17 at 16:14
  • You have still mixed up out_layer and dst_layer. – AndreJ Mar 11 '17 at 18:00
  • 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... – dgj32784 Mar 11 '17 at 21:41
  • I don't have much experience with python, but do you tell the driver to create a polygon layer anywhere? – AndreJ Mar 12 '17 at 6:32
  • 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) – dgj32784 Mar 12 '17 at 16:30

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.