If you are a keen GDAL user, you have probably noticed by now that you can extract a rectangular portion of an image into a new image. Sometimes you want to go a little further and actually extract only the cells within a polygonal area from a shapefile or similar vector file. In this article I will show you how.

Overview

There are basically four steps to follow:

- create a shpfile containing only the features you wish to use as your clip mask
- determine the extent of your clipmask
- perform a rectangular clip on the image
- perform a polygonal clip on the image

Heres a little preview of what we will be able to produce:

A quick diversion to add support for MrSid

If my previous article on mosaicking I illustrated how GDAL can can be built with ECW support. Before I get started with my clipping, I am going to add MrSid support since some of the rasters I plan to clip are in MrSid format. To do this you need to download the MrSid geo express sdk from their website (requires registration). Once you have it downloaded, unpack it into /usr/local/GeoExpressSDK and then follow the gdal build guide from my mosaicking article while adding a couple of options to the configure line:

./configure --with-libtiff=internal --with-geotiff=internal \
--with-ecw=/usr/local --with-mrsid=/usr/local/GeoExpressSDK --without-jp2mrsid
make
sudo make install

Once installed, you can check if your gdal supports MrSid by doing:

/usr/local/bin/gdalinfo --formats | grep SID

Which should result in a line like this being displayed:

  MrSID (ro): Multi-resolution Seamless Image Database (MrSID)

Create your mask shapefile

For my purposes, I want to clip data out for the Eastern Cape Province of South Africa. So I loaded a national dataset of SA Provinces in QGIS, selected the Eastern Cape province and the exported that to its own shapefile by right clicking on its legend entry and choosing the ‘Save selection as shapefile’ option. It is important that your shapefile contain only the features needed for clipping as we will be using its bounding box as part of the clip process.

Determine the extent of your clip mask

With a few lines of bash and the ogrinfo programme we can determine the bounding box of a shapefile e.g.:

#!/bin/bash
SHPFILE=$1
BASE=`basename $SHPFILE .shp`
EXTENT=`ogrinfo -so $SHPFILE $BASE | grep Extent \
| sed 's/Extent: //g' | sed 's/(//g' | sed 's/)//g' \
| sed 's/ - /, /g'`
EXTENT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
echo $EXTENT

If we run the script it produces output like this:

22.735740 -30.001800 30.194720 -34.213860

This gives you the extents in the order and format expected by gdal_translate which we need for our clip operation.

Clip the raster to the bounding box of your shapefile

To do the clip to bounding box, we can simply add a few more lines to our script:

RASTERFILE=$2
gdal_translate -projwin $EXTENT -of GTiff $RASTERFILE /tmp/`basename $RASTERFILE .sid`_bbclip.tif

Ok so that should have clipped out your raster to the bounding box of your shapefile.

Perform a polygonal clip on the image

Last on our list is to clip out just the polygonal area from the clipped-to-bounding-box version of our image. To do this we use gdalwarp e.g.:

gdalwarp -co COMPRESS=DEFLATE -co TILED=YES -of GTiff \
-r lanczos -cutline $SHPFILE \
/tmp/`basename $RASTERFILE .sid`_bbclip.tif /tmp/`basename $RASTERFILE .sid`_shpclip.tif

Tying it all together

Lets see what our final script looks like then:

#!/bin/bash
SHPFILE=$1
BASE=`basename $SHPFILE .shp`
EXTENT=`ogrinfo -so $SHPFILE $BASE | grep Extent \
| sed 's/Extent: //g' | sed 's/(//g' | sed 's/)//g' \
| sed 's/ - /, /g'`
EXTENT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
echo "Clipping to $EXTENT"
RASTERFILE=$2
gdal_translate -projwin $EXTENT -of GTiff $RASTERFILE /tmp/`basename $RASTERFILE .sid`_bbclip.tif
gdalwarp -co COMPRESS=DEFLATE -co TILED=YES -of GTiff \
-r lanczos -cutline $SHPFILE \
/tmp/`basename $RASTERFILE .sid`_bbclip.tif /tmp/`basename $RASTERFILE .sid`_shpclip.tif

Save the script and make it executable. I saved mine to mine as ~/bin/clip.sh and then made it executable using

chmod +x ~/bin/clip.sh

To invoke the script you would do something like:

~/bin/clip myshapefile.shp mybigraster.sid

Of course there are some things that can be improved in the script – for example its hard coded to assume that you are working with MrSid data. However the script should provide you with a starting point to clipping out your own raster data using arbitary polygons.