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.
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.
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
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.
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.
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