Integrate the Python Script that Clips Multiple Rasters at Once in the QGIS Processing Framework (SEXTANTE)


This post is about how to integrate a Python script for clipping rasters in the SEXTANTE toolbox for QGIS (called Processing framework in QGIS 2.0). The script is not run inside the Python Console, but it has a specific interface, and it can be used in the Modeler in order to create automated workflows.
This script allows for the clipping of multiple rasters using a shapefile boundary or the coordinates of a rectangle (similarly to this post). Clipping is useful before processing rasters, in order to limit the classification to the study area.
A video tutorial is at the end of this post.

Interface of the script Clip Multiple Rasters
It is possible to create Python scripts that integrate in the SEXTANTE toolbox using the function in the SEXTANTE toolbox > Scripts > Tools > Create new script, and writing the code. Therefore, it is possible to create workflows for GIS analyses with a customized interface for input and output data. For more information about the parameters see here.

The following is the code of this script:
##Raster=group
##Input=multiple raster
# upper left point
##UX=number 0
##UY=number 0
# lower right point
##LX=number 0
##LY=number 0
##Use_shapefile=boolean
##Shapefile=vector
##No_data=boolean
##No_data_value=number 0
##Output_directory=folder

import os
import sextante

lddLrs = Input.split(';')

if Use_shapefile == 0:
if No_data == 0:
for lyr in lddLrs:
os.system("gdal_translate -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))
else:
for lyr in lddLrs:
os.system("gdal_translate -a_nodata " + str(No_data_value) + " -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))
else:
if No_data == 0:
for lyr in lddLrs:
os.system("gdalwarp -cutline " + str(Shapefile) + " -crop_to_cutline -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))
else:
for lyr in lddLrs:
os.system("gdalwarp -dstnodata " + str(No_data_value) + " -cutline " + str(Shapefile) + " -crop_to_cutline -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))

How this script works:
The first part of the script sets the parameters that are shown in the interface (the double # are required):
##Raster=group
##Input=multiple raster
# upper left point
##UX=number 0
##UY=number 0
# lower right point
##LX=number 0
##LY=number 0
##Use_shapefile=boolean
##Shapefile=vector
##No_data=boolean
##No_data_value=number 0
##Output_directory=folder

Where:
  • Input are the selected rasters.
  • UX is the X coordinate of the upper left point that defines the rectangle.
  • UY is the Y coordinate of the upper left point that defines the rectangle.
  • LX is the X coordinate of the lower left point that defines the rectangle.
  • LY is the Y coordinate of the lower left point that defines the rectangle.
  • Use shapefile is an option. If true, then the Shapefile is used for clipping.
  • Shapefile is the boundary used for clipping the rasters.
  • No data is an option. If true, then the No data value is used.
  • No data value is a raster value Pixels with this value are set to No data.
  • Output directory is where the clipped rasters are saved. Clipped rasters have the name clip_ + original raster name and are saved as .tif files.

The second part is the actual script, that uses gdal_translate for clipping by coordinates, or gdalwarp for clipping by shapefile.
import os
import sextante

lddLrs = Input.split(';')

if Use_shapefile == 0:
if No_data == 0:
for lyr in lddLrs:
os.system("gdal_translate -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))
else:
for lyr in lddLrs:
os.system("gdal_translate -a_nodata " + str(No_data_value) + " -projwin " + str(UX) + " " + str(UY) + " " + str(LX) + " " + str(LY) + " -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))
else:
if No_data == 0:
for lyr in lddLrs:
os.system("gdalwarp -cutline " + str(Shapefile) + " -crop_to_cutline -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))
else:
for lyr in lddLrs:
os.system("gdalwarp -dstnodata " + str(No_data_value) + " -cutline " + str(Shapefile) + " -crop_to_cutline -of GTiff " + str(lyr) + " " + Output_directory + "/clip_" + str(os.path.basename(str(lyr))))
sextante.load(str(Output_directory + "/clip_" + str(os.path.basename(str(lyr)))))

You can download the script and the help file for QGIS 1.8 from here (or for QGIS 2 from here).

How to use this script:
The following example uses a sample dataset (a Landsat 8 image available from the U.S. Geological Survey) that you can download here and a shapefile from here.

  1. Copy the script files in the directory scripts
    • Unzip the downloaded file and copy the two files inside the script directory in your home/user directory (that is .qgis/sextante/scripts for QGIS 1.8 or .qgis2/processing/scripts for QGIS 2);
  1. Load the rasters and the shapefile in QGIS
    • Load the bands you wish to clip (e.g. 2, 3, and 4) and the mask shapefile in QGIS;
  1. Clip the rasters
    • In the SEXTANTE toolbox double click on Scripts > Raster > Clip Multiple Rasters; select the Input rasters; under Use shapefile select Yes (or True); select the shapefile (e.g. mask); optionally, under No data select Yes and set the No data value (e.g. to 0); select an output directory and click OK; 
    • After a while, the clipped rasters will be loaded in QGIS.

Attention! This script for QGIS 2 does not load the clipped rasters in QGIS automatically.

Clipping by coordinates works very similarly; you need to type the coordinates that define the rectangle under  UX, UY, LX, LY. Unfortunately, SEXTANTE for QGIS 1.8 has an issue that precludes the typing of long coordinates (solved in QGIS 2), therefore you should use the script that I have described here, or use the script in QGIS 2.
It is worth mentioning that a similar tool is implemented in the development version of QGIS 2.1, in the Processing toolbox > GDAL > ClipByExtent (see here).

The following is the video of this tutorial.

Newer posts Older posts