Here's a command line method using gdalwarp, 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}"
gdalwarp