5

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 ?

  • The easiest way is to use a scripted solution, where you first split your polygon into features and then crop your raster element-wise. What operating system are you using? – Curlew Mar 31 '13 at 16:39
  • I try to develop a multi-platform QGIS plugin so I use MacOS, Windows and Linux (Ubuntu). – Lambdacrash Mar 31 '13 at 16:41
  • If you can write scripts then what is the problem? Do you want someone to write a code for you? Aren't you able to clip the same raster 20 times in the cycle? – SS_Rebelious Mar 31 '13 at 16:45
  • I don't want to clip my raster 20 times. I have a shapefile containing x polygons. I want to clip my initial raster by the first polygon and get a new raster output, then, clip my initial raster by the second polygon and get new raster output, and so on ... I can't find how to do this with gdalwarp – Lambdacrash Mar 31 '13 at 16:52
  • @Lambdacrash, pseudocode: (for polygon in layer: clip_initial_raster.by_polygon, save output). – SS_Rebelious Mar 31 '13 at 17:03
5

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}"

0

I found the first part of the solution here. I just have to iterate over features.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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