Mask raster by extent in Python using rasterio

Question

I want to clip one raster based on the extent of another (smaller) raster. First I determine the coordinates of the corners of the smaller raster using

import rasterio as rio
import gdal
from shapely.geometry import Polygon

src = gdal.Open(smaller_file.tif)
ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()
lrx = ulx + (src.RasterXSize * xres)
lry = uly + (src.RasterYSize * yres)
geometry = [[ulx,lry], [ulx,uly], [lrx,uly], [lrx,lry]]

This gives me the following output geometry = [[-174740.0, 592900.0], [-174740.0, 2112760.0], [900180.0, 2112760.0], [900180.0, 592900.0]]. (Note that the crs is EPSG: 32651). Now I would like to clip the larger file using rio.mask.mask(). According to the documentation, the shape variable should be GeoJSON-like dict or an object that implements the Python geo interface protocol (such as a Shapely Polygon). Therefore I create a Shapely Polygon out of the variable geometry, using

roi = Polygon(geometry)

Now everything is ready to use the rio.mask() function.

output = rio.mask.mask(larger_file.tif, roi, crop = True)

But this gives me the following error

TypeError: 'Polygon' object is not iterable

What do I do wrong? Or if someone knows a more elegant way to do it, please let me know.

(Unfortunately I cannot upload the two files since they're too large)

Answer 1

I found your question when I needed to figure out this kind of clipping myself. I got the same error and fixed it the following way:

rasterio.mask expects a list of features, not a single geometry. So the algorithm wants to run masking over several features bundled in an iterable (e.g. list or tuple) so we need to pass it our polygon within a list (or tuple) object. The code you posted works after following change:

roi = [Polygon(geometry)]

All we have to do is to enclose the geometry in a list/tuple and then rasterio.mask works as expected.