2013/09/11

How to Clip Multiple Rasters at Once Using QGIS, Python and GDAL

This is a tutorial about how to clip multiple rasters at once in QGIS, running a very simple Python script that uses the command gdal_translate.
This way, every raster loaded in QGIS will be clipped to the same area, defined by a rectangle.
A video tutorial is also available at the end of this post.

The Python script is the following:

# input variables
# where to save clipped rasters
outputDir = "/home/user/Desktop/clip"
# upper left point
UX, UY =  298545 , 4628145
# lower right point
LX,  LY = 308715 , 4622745

# actual script
import os
lddLrs = qgis.utils.iface.legendInterface().layers()
for lyr in lddLrs:
if (lyr.type()==QgsMapLayer.RasterLayer):
os.system("gdal_translate -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr.source()) + " " + outputDir + "/" + str(lyr.name()) + ".tif")
qgis.utils.iface.addRasterLayer(str(outputDir + "/" + lyr.name() + ".tif"), str(lyr.name() + "_clip.tif"))

How this script works:
# where to save clipped rasters
outputDir = "/home/user/Desktop/clip"
# upper left point
UX, UY =  298545 , 4628145
# lower right point
LX,  LY = 308715 , 4622745
The input variables are: the directory where to save clipped rasters, and the XY coordinates that define the clip rectangle: U is the upper left point and L the lower right point. Coordinates can be copied using ctrl + c from the QGIS current map coordinates or from the Identify Results.

import os
Import the required module for using gdal_translate

lddLrs = qgis.utils.iface.legendInterface().layers()
Get all the layers loaded in QGIS

for lyr in lddLrs:
Loop through all the loaded layers

if (lyr.type()==QgsMapLayer.RasterLayer):
Check if the layer is a raster

os.system("gdal_translate -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr.source()) + " " + outputDir + "/" + str(lyr.name()) + ".tif")
qgis.utils.iface.addRasterLayer(str(outputDir + "/" + lyr.name() + ".tif"), str(lyr.name() + "_clip.tif"))
For each raster, use gdal_translate to clip the raster using the input coordinates, and save it as .tif file in the output directory

Hint: it is possible to add a no data value (e.g. 0) to the clipped raster, by adding the code
 -a_nodata 0 
after gdal_translate

qgis.utils.iface.addRasterLayer(str(outputDir + "/" + lyr.name() + ".tif"), str(lyr.name() + ".tif"))
Load the clipped rasters in QGIS


How to use this script:
  1. Load the rasters to be clipped in QGIS;
  2. Copy this script in a text editor; change the input variables according to your needs; select all and copy the text;
  3. In QGIS, open the Plugins > Python Console; paste the code in the Python Console and press Enter twice.

Probably, I am going to add a similar script to the Semi-Automatic Classification Plugin as a new utility for clipping multiple rasters.
The following is the video of this tutorial using a sample dataset (a Landsat 8 image available from the U.S. Geological Survey) that you can download here.

No comments:

Post a Comment