Polygonizing global raster using Python/GDAL?

Question

I'm looking to polygonize a global raster with 29 different categories. However, when I try to run it in QGIS, it either crashes or takes forever to load.

I'm now trying to run it in python.

import subprocess, glob
script = "/usr/bin/gdal_polygonize.py"
for raster in glob.glob("/home/chiara/Desktop/WorldSoils/*.tif"):
 outshp = raster[:-4] + ".shp"
 subprocess.call(["python", script, raster, "-f", "ESRI Shapefile", outshp])

However, this still takes forever to load. Is there a way to polygonize on this raster by splitting up the earth into several grid cells (perhaps each about 6-degrees) and polygonizing each of those grid cells one at a time?

Answer 1

You can indeed split a raster into smaller tiles, but you will need to know the extents of each tile. gdal_translate and -projwin will give you the result you want.

Another solution, one I have used on large rasters, is to simplify/re-clasify the raster reducing its size. I have done this with gdal_calc, outputing a GeoTiff, --type Byte, with the following creation options --co "NUM_THREADS=ALL_CPUS" --co "TILED=YES" --co "COMPRESS=LZW" --co "SPARSE_OK=TRUE", --NoDataValue 00.000 and by simplifying the values to a range of values, in your case, 1 - 29. I have found that polygonizing these can be a lot quicker.