Python GDAL: Clip rasterfile with multi-polygon shapefile

Question

I am trying to clip a raster file with a multi-polygon shapefile to get separate clipped-raster file for each polygon with its 'id' as part of the name of clipped-raster file. I tried the following code by it throws an error.

code:

import ogr
import subprocess

shapefile = ("/shape.shp")
in_raster = ("/whole.tif")

ds = ogr.Open(shapefile)
lyr = ds.GetLayer(0)

lyr.ResetReading()
ft = lyr.GetNextFeature()

while ft:

   shape_id= ft.GetFieldAsString('id')

    out_raster = in_raster.replace('.tif', '_%s.tif' % shape_id.replace(' ', '_'))    
    subprocess.call(['gdalwarp', in_raster, out_raster, '-cutline', shapefile, 
                     '-crop_to_cutline', '-cwhere', "'id'='%s'" % shape_id])

    ft = lyr.GetNextFeature()
ds = None

Error:
ERROR 1: Did not get any cutline features

Answer 1

It's often easier to use the GDAL python wrappers. See my attempt below. If you still get errors, check that both raster and shapefile have geographic information (i.e. you can overlap them in QGIS or something like that).

from osgeo import gdal
from osgeo import ogr

def dump_poly(raster_fname, vector_fname, ifeat):
    # New filename. Assumes input raster file has '.tif' extension
    # Might need to change how you build the output filename
    fname_out = raster_fname.replace(".tif", f"_{int(ifeat):05}.tif")
    # Do the actual clipping
    g = gdal.Warp(fname_out, raster_fname, format="GTiff",
                    cutlineDSName=vector_fname,
                    cutlineWhere=f"id={ifeat:d}",
                    cropToCutline=True)
    # Return the fname just because
    return fname_out


vector_fname = "vector_file.shp"
raster_fname = "raster_file.tif"
# Open raster file, select first and unique layer
gg = ogr.Open(vector_fname)
layer = gg.GetLayerByIndex(0)
# Loop over all features
for ifeat in layer:
    ifeat = GetFieldAsInteger('id')
    fname_out = dump_poly(raster_fname, vector_fname, ifeat)
    print(fname_out)