[GIS] How to clip raster by multiple polygons in multi rasters

gdalpythonqgisraster

Using QGIS, I want to clip a raster by multiple polygons contained in a single shapefile. I read that gdalwarp can clip my raster by multiple polygons BUT the output is a single merged clipped raster.

I want to have one clipped raster per polygon. If I have 20 polygons, I want to have 20 clipped rasters at the end 🙂

I can't find any option of gdalwarp that allows me to do this. Do you have any idea ?

Best Answer

Here's a command line method using , formulated as a Windows batch file:

@echo off
setlocal
set index=D:\source\nts_index_250k.shp
set tiles=117D 117A 116O 116P

for %%a in (%tiles%) do (
  @echo gdalwarp -cutline %index% ^
  -cwhere "'TILE_NAME' = '%%a'" ^
  -cblend 3 ^
  %1 %2_%%a.tif
  )
endlocal
goto :eof

Usage

crop_by_poly.bat infile.tif

Explanation

As written, the batchfile does no work, just prints to console the commands needed to do the work (remove @echo to make it execute instead of print):

gdalwarp -cutline D:\source\nts_index_250k.shp   -cwhere "'TILE_NAME' = '117D'"   -cblend 3 dem_20m_shade.tif 117D_dem_20m_shade.tif
gdalwarp -cutline D:\source\nts_index_250k.shp   -cwhere "'TILE_NAME' = '117A'"   -cblend 3 dem_20m_shade.tif 117A_dem_20m_shade.tif

^: line continuation character

index: polygon data source, in the case Canadian National Topographic System tiles

tiles: values to select from the polygon feature records. Attribute table looks like:

+-----------+-------------------+
| TILE_NAME |       NAME        |
+-----------+-------------------+
| 117C      | Demarcation Point |
| 117D      | Herschel Island   |
| 107C      | Mackenzie Delta   |
+-----------+-------------------+

-cwhere "'TILE_NAME' = '%%a'" ^ : select this or these features from source, becomes -cwhere "'TILE_NAME' = 117D'" etc. Be mindful of quotes, outermost must be double and inner must be single, and always in pairs.

-cblend 3: optional, expand crop area by 3 pixels, to allow for future seamless overlap/mosaicking

%1 %2: becomes infile.tif outprefix_infile.tif, so in this example 117D_infile.tif, 117A_infile.tif, ...

The output image will have same extent as original raster with the clipped out regions marked as nodata. Add the -crop_to_cutline option to have set the extent to match the clipping polygon, but be aware that the edge matching is not perfect at this time when used with cblend (gdal issue #5981)

A python solution was asked for, the logic illustrated with this batchfile could be adapted to python by incorporating For looping folder to batch clip rasters by polygon using python and QGIS?

Untested and not elegant or pythonic, but should get someone started:

from subprocess import call
clipper = r'D:\source\nts_index_250k.shp'
tiles = '117D 117A 116O 116P'.split()
for tile in tiles:
    warp = '''gdalwarp -cutline {clipper} -cwhere "'TILE_NAME' = '{tile}'" {infile} {outfile}'''.format(clipper=clipper,
        tile=tile,
        infile=infile,
        outfile=tile + infile
        )
    call(warp)

Dear lazy web: The ideal answer to this question would be to operate from inside qgis: select a polygon layer (or selection polys within layer), select a raster, and run "clip_for_each {selected} {outfolder}"

Related Question