Thursday 28 March 2019

Using GDAL polygonize from Python?


I have been using gdal from the command line to convert an asc file to a GeoJSON output. I can do this successfully:


gdal_polygonize.py input.asc -f "GeoJSON" output.json


Now I wish to use Python and follow this process for a range of files.


import gdal
import glob
for file in glob.glob("dir/*.asc"):

new_name = file[:-4] + ".json"
gdal.Polygonize(file, "-f", "GeoJSON", new_name)

Hpwever, for exactly the same file I get the following error TypeError: in method 'Polygonize', argument 1 of type 'GDALRasterBandShadow *'


Why does the command line version work and the python version not?



Answer



You are confusing the use of the gdal_polygonize command line utility with the python function gdal.Polygonize().


As you mentioned, you've managed to use the command line utility successfully; however, the Python function works differently and expects different arguments than those specified in the utility. The first argument should be a GDAL Band object, not a string, so this is why you get your error.


To get the Band object you need to open the input file using gdal.Open() and use the GetRasterBand() method to get your intended band. Additionally, you need to create an output layer in which the resulting polygons will be created.


The Python GDAL/OGR Cookbook has a good example on how to use this function. The required parameters are explained in a bit more detail here.





Alternative


Following on from your comment, if you would prefer to keep using the command line utility, one solution is to call it from within a Python script using the subprocess module i.e.


import subprocess

script = 'PATH_TO_GDAL_POLYGONIZE'
for in_file in glob.glob("dir/*.asc"):
out_file = in_file[:-4] + ".json"
subprocess.call(["python",script,in_file,'-f','GeoJSON',out_file])


This loops through the files and updates the input/output paths.


This way you get the result that you are use to getting from the utility, but the ability to loop through your files. If you do plan on going this way, the subprocess documention will be useful.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...