0

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
3
  • Print and add to your question what you send to subprocess.call. It should be something that you can run successfully with command line gdalwarp.
    – user30184
    Jan 21, 2022 at 7:08
  • It might just be the path to the shapefile that is wrong. If the shapefile is in the same directory as your python script you can try with shapefile = 'shape.shp'
    – Dataform
    Jan 21, 2022 at 7:56
  • @Dataform, I tried but nothing happened, same result Jan 21, 2022 at 11:40

1 Answer 1

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)
2
  • worked with some little modification, as : ifeat = GetFieldAsString('id') to ifeat = GetFieldAsInteger('id'). Because my id field has integers. Thanks a lot. I was struggling with it for last 3 days. Jan 21, 2022 at 17:07
  • 1
    Ooops! Well spotted, have now updated the code with your fix. Glad it worked.
    – Jose
    Jan 24, 2022 at 11:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.