| |

How to Clip a Raster to a Polygon or an Extent with Python (Extract by Mask)

Clipping, extracting to masks, and subsetting rasters are common GIS operations. Using the gdal package for Python it is very easy to automate these processes. This tutorial will demonstrate how to use the Warp() function from the gdal Python API to clip rasters to a specified extent and to clip rasters to a polygon layer. These operations are also often referred to as mask extraction or extract to mask.

Let’s get started.

Install GDAL

You’ll need to have gdal installed in your Python environment to complete this tutorial. Unfortunately, gdal is not the easiest package to install, especially if you are using a non-anaconda Python distribution. But don’t worry. I’ve written a guide to installing gdal that you can follow. I do recommend using an Anaconda Python distribution because it makes gdal installation and updates much easier. I prefer to use miniconda because it is lightweight. You can follow this guide for how to install gdal with anaconda/miniconda.

Input Data

Now we need some input data to work with. I’m going to be using a 1 m lidar-derived DEM from the USGS. I’ll also be using a polygon I created that overlaps a portion of the raster image. You can download and create these data yourself, or you can sign up for my full GDAL course which provides step-by-step video instruction, code, and data. Once you have the data, you’re ready to get started.

Clip Raster to an Extent

Start by importing gdal and osr (automatically installed with gdal. Depending on how you installed gdal you might have a slightly different import statement than the one I am using.

Specify the input raster file and output raster file. Then specify the bounds to clip the raster to. Make sure these bounds are within the bounds of the input raster. If they are not you will receive an error. Also, be sure your bounds are in the same spatial reference system (SRS) as your input raster. It is possible to clip to bounds with a different SRS, but we will not cover that here. You can learn more in the full course.

from osgeo import gdal, osr

fn_in = r"path/to/input.tif"
fn_extent = 'path/to/output.tif'

bounds_utm11 = (641000., 4851000., 649000., 4859000.)

When performing clips and extracts I like to specify the output coordinate system, just to be safe. Here the coordinate system is NAD83 UTM Zone 11 N. I use osr to set the spatial reference below.

utm11_srs = osr.SpatialReference()
utm11_srs.SetWellKnownGeogCS("NAD83")
utm11_srs.SetUTM(11, 1)

Now we’re ready to perform the actual clip ( or extraction) with gdal.Warp(). Specify the output file name, the input file name, the output bounds, and the SRS of the output bounds. This is done with the following line of Python code.

gdal.Warp(fn_extent, fn_in, outputBounds=bounds_utm11, outputBoundsSRS=utm11_srs)

Run the full script and open the output raster to see the result. If you are using the data from my course, you will get a result like the one below.

Original raster
Clipped to extent

Clip Raster to Polygon

We will also use gdal.Warp() to clip a raster to a polygon with Python. To do this we’ll need to set things up just a little bit differently.

To start, import gdal. There is no need to import osr for this example because we’re not specifying coordinates. Also, set up the input raster file name, the polygon file name, and the output raster file name. I’ll be using the same input raster as the example above and clipping with the polygon shown in the image below. The polygon can be in any file format that is readable by ogr. Here I’m using a geopackage.

from osgeo import gdal

fn_in = r"path/to/input.tif"
fn_clip = 'path/to/output.tif'
fn_poly = r"path/to/polygon.gpkg"

Now we’ll pass the appropriate arguments to gdal.Warp() to clip the raster to a polygon. We need to specify the input file name, output file name, polygon dataset file name, polygon dataset layer, and specify to crop to the polygon. That can be done in a single line of code. Note that ‘polygon’ is the name of the layer inside the geopackage that I’m clipping with.

gdal.Warp(fn_clip, fn_in, cutlineDSName=fn_poly, cutlineLayer='polygon', cropToCutline=True)

After running the code your raster will be clipped to the polygon shape. You can see my result in the images below. You can also specify a where clause in gdal.Warp() to select specific features from a polygon layer. For example, I could use that clause to select only one of the polygons in the layer below instead of clipping to both. A full demonstration of this is available in the full GDAL course.

Input raster and polygon
Clipped to polygon

Do more with GDAL

Hopefully, this tutorial has helped you learn one way that gdal can help you automate GIS processes with Python. Coding up GIS tasks makes them exactly reproducible and can save you hours of time. With GDAL you can also resize rasters, merge (or mosaic) multiple rasters together, perform terrain analysis and interpolation, and do many other raster operations. If you’re interested in learning more about how to improve your geospatial skills check out the full GDAL course that is linked in this article, and the other course offerings from geospatialschool.com by following the link below.

Whether you’re looking to take your GIS skills to the next level, or just getting started with GIS, we have a course for you! We’re constantly creating and curating more courses to help you improve your geospatial skills.

All of our courses are taught by industry professionals and include step-by-step video instruction so you don’t get lost in YouTube videos and blog posts, downloadable data so you can reproduce everything the instructor does, and code you can copy so you can avoid repetitive typing

Similar Posts