Clipping rasters can be trivial with a desktop GIS like QGIS or with command line tools like GDAL. However, I recently ran into a situation where I needed to clip large rasters in an automated, online Python process. It simply wouldn’t do to interrupt the procedure and clip them myself. The Python bindings for GDAL/OGR are pretty neat but they are very low-level; how could I use Python to clip the rasters as part of a continuous Python session?

Jared Erickson posted an excellent tutorial on this topic, one of many in his “Python GDAL/OGR Cookbook”. Even this simple example hints at how complicated something like clipping can be with the low-level GDAL/OGR API. However, there are still some limitations to the example he provided:

  • It only works for single clipping features.
  • It only works for multiple-band rasters.
  • The clipping features must be bounded within the raster.

Here I’ll present a solution to the second and third points: Clipping a raster, stacked or single-band, with clipping features that extend beyond its bounds.

Improvement: Support for Single-Band Images

This is an easy fix. We simply need to use a try...catch sequence of statements to anticipate the error in our NumPy indexing when a single-band raster is presented.

# Multi-band image?
try:
    clip = rast[:, ulY:lrY, ulX:lrX]

# Nope: Must be single-band
except IndexError:
    clip = rast[ulY:lrY, ulX:lrX]

Improvement: Out-of-Bounds Clipping Features

The second part is not so easy. The clipping in Jared’s example is based on NumPy arrays and therefore must be conformal. One thing that can happen when the clipping features extend “above” a raster’s extent is that the estimated pixel coordinate for the upper-left corner become negative—and negative pixels don’t make sense! We have to catch this case, remember the negative pixel coordinate, and set it to zero.

# If the clipping features extend out-of-bounds and ABOVE the raster...
if gt[3] < maxY:
    # In such a case... ulY ends up being negative--can't have that!
    iY = ulY
    ulY = 0

However, in doing this, the clipping features are effectively “pushed down” (southwards) because we prevented the upper-left corner coordinate from being negative. To compensate, we must pull the clipping features “back up” (northwards). Note that here we need to add an else clause for the case when the clipping features don’t extend above the raster.

# If the clipping features extend out-of-bounds and ABOVE the raster...
if gt[3] < maxY:
    # The clip features were "pushed down" to match the bounds of the
    #   raster; this step "pulls" them back up
    premask = image_to_array(raster_poly)
    # We slice out the piece of our clip features that are "off the map"
    mask = np.ndarray((premask.shape[-2] - abs(iY), premask.shape[-1]), premask.dtype)
    mask[:] = premask[abs(iY):, :]
    mask.resize(premask.shape) # Then fill in from the bottom

    # Most importantly, push the clipped piece down
    gt2[3] = maxY - (maxY - gt[3])

else:
    mask = image_to_array(raster_poly)

Finally, we have to deal with the converse problem: clipping features that extend below the raster’s bounds. This is trickier. What I’ve done is created a larger NumPy array onto which the clipping features are projected.

# Clip the image using the mask
try:
    clip = gdalnumeric.choose(mask, (clip, nodata))

# If the clipping features extend out-of-bounds and BELOW the raster...
except ValueError:
    # We have to cut the clipping features to the raster!
    rshp = list(mask.shape)
    if mask.shape[-2] != clip.shape[-2]:
        rshp[0] = clip.shape[-2]

    if mask.shape[-1] != clip.shape[-1]:
        rshp[1] = clip.shape[-1]

    mask.resize(*rshp, refcheck=False)

    clip = gdalnumeric.choose(mask, (clip, nodata))

The Complete Source Code

Here’s a complete Python function for clipping a raster. It follows Jared’s example but expands on the comments throughout and includes the two improvements I’ve described. It does not yet support clipping by multiple features at once. Also, I haven’t evaluated (nor come across the problem of) its performance with clipping features that extend only to the side of the raster; not sure if this would be a problem.

The function’s API consists of two required arguments: The raster to clip and the file path to the clipping features (e.g., a Shapefile). An optional GDAL GeoTransform can be provided. Also, the NoData value can be set.

from osgeo import gdal, gdalnumeric, ogr
from PIL import Image, ImageDraw
import os
import numpy as np

def clip_raster(rast, features_path, gt=None, nodata=-9999):
    '''
    Clips a raster (given as either a gdal.Dataset or as a numpy.array
    instance) to a polygon layer provided by a Shapefile (or other vector
    layer). If a numpy.array is given, a "GeoTransform" must be provided
    (via dataset.GetGeoTransform() in GDAL). Returns an array. Clip features
    must be a dissolved, single-part geometry (not multi-part). Modified from:

    http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html
    #clip-a-geotiff-with-shapefile

    Arguments:
        rast            A gdal.Dataset or a NumPy array
        features_path   The path to the clipping features
        gt              An optional GDAL GeoTransform to use instead
        nodata          The NoData value; defaults to -9999.
    '''
    def array_to_image(a):
        '''
        Converts a gdalnumeric array to a Python Imaging Library (PIL) Image.
        '''
        i = Image.fromstring('L',(a.shape[1], a.shape[0]),
            (a.astype('b')).tostring())
        return i

    def image_to_array(i):
        '''
        Converts a Python Imaging Library (PIL) array to a gdalnumeric image.
        '''
        a = gdalnumeric.fromstring(i.tobytes(), 'b')
        a.shape = i.im.size[1], i.im.size[0]
        return a

    def world_to_pixel(geo_matrix, x, y):
        '''
        Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
        the pixel location of a geospatial coordinate; from:
        http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#clip-a-geotiff-with-shapefile
        '''
        ulX = geo_matrix[0]
        ulY = geo_matrix[3]
        xDist = geo_matrix[1]
        yDist = geo_matrix[5]
        rtnX = geo_matrix[2]
        rtnY = geo_matrix[4]
        pixel = int((x - ulX) / xDist)
        line = int((ulY - y) / xDist)
        return (pixel, line)

    # Can accept either a gdal.Dataset or numpy.array instance
    if not isinstance(rast, np.ndarray):
        gt = rast.GetGeoTransform()
        rast = rast.ReadAsArray()

    # Create an OGR layer from a boundary shapefile
    features = ogr.Open(features_path)
    if features.GetDriver().GetName() == 'ESRI Shapefile':
        lyr = features.GetLayer(os.path.split(os.path.splitext(features_path)[0])[1])

    else:
        lyr = features.GetLayer()

    # Get the first feature
    poly = lyr.GetNextFeature()

    # Convert the layer extent to image pixel coordinates
    minX, maxX, minY, maxY = lyr.GetExtent()
    ulX, ulY = world_to_pixel(gt, minX, maxY)
    lrX, lrY = world_to_pixel(gt, maxX, minY)

    # Calculate the pixel size of the new image
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    # If the clipping features extend out-of-bounds and ABOVE the raster...
    if gt[3] < maxY:
        # In such a case... ulY ends up being negative--can't have that!
        iY = ulY
        ulY = 0

    # Multi-band image?
    try:
        clip = rast[:, ulY:lrY, ulX:lrX]

    except IndexError:
        clip = rast[ulY:lrY, ulX:lrX]

    # Create a new geomatrix for the image
    gt2 = list(gt)
    gt2[0] = minX
    gt2[3] = maxY

    # Map points to pixels for drawing the boundary on a blank 8-bit,
    #   black and white, mask image.
    points = []
    pixels = []
    geom = poly.GetGeometryRef()
    pts = geom.GetGeometryRef(0)

    for p in range(pts.GetPointCount()):
        points.append((pts.GetX(p), pts.GetY(p)))

    for p in points:
        pixels.append(world_to_pixel(gt2, p[0], p[1]))

    raster_poly = Image.new('L', (pxWidth, pxHeight), 1)
    rasterize = ImageDraw.Draw(raster_poly)
    rasterize.polygon(pixels, 0) # Fill with zeroes

    # If the clipping features extend out-of-bounds and ABOVE the raster...
    if gt[3] < maxY:
        # The clip features were "pushed down" to match the bounds of the
        #   raster; this step "pulls" them back up
        premask = image_to_array(raster_poly)
        # We slice out the piece of our clip features that are "off the map"
        mask = np.ndarray((premask.shape[-2] - abs(iY), premask.shape[-1]), premask.dtype)
        mask[:] = premask[abs(iY):, :]
        mask.resize(premask.shape) # Then fill in from the bottom

        # Most importantly, push the clipped piece down
        gt2[3] = maxY - (maxY - gt[3])

    else:
        mask = image_to_array(raster_poly)

    # Clip the image using the mask
    try:
        clip = gdalnumeric.choose(mask, (clip, nodata))

    # If the clipping features extend out-of-bounds and BELOW the raster...
    except ValueError:
        # We have to cut the clipping features to the raster!
        rshp = list(mask.shape)
        if mask.shape[-2] != clip.shape[-2]:
            rshp[0] = clip.shape[-2]

        if mask.shape[-1] != clip.shape[-1]:
            rshp[1] = clip.shape[-1]

        mask.resize(*rshp, refcheck=False)

        clip = gdalnumeric.choose(mask, (clip, nodata))

    return (clip, ulX, ulY, gt2)