Precicely cropping a raster in rasterio

You can use gdalwarp or raster.merge to crop a raster to specified bounds, optionally reprojecting and modifying the resolution in the process.

But the bounds of the output may differ slightly from the requested bounds due to rounding.

You can crop (or expand) a raster to exactly the desired bounds in rasterio, it just takes a bit more work.

Take this raster (which by default will be in UTM projection, buffered slightly larger than the requested bounds)

bounds = rasterio.coords.BoundingBox(-121.6271, 37.1622, -121.6184, 37.1650)
url = f"https://api.gpxz.io/v1/elevation/hires-raster?key={GPXZ_API_KEY}&bbox_left={bounds.left}&"
url += f"bbox_right={bounds.right}&bbox_bottom={bounds.bottom}&bbox_top={bounds.top}"

You’d like to crop it to the exact latlon bounds, and also reproject it to WGS84 latlon.

dst_crs = "EPSG:4326"

Rasterio’s reprojection function can warp a raster to a projection, bounds, and resolution, as long as you have an affine transformation and a desired shape.

To get the transform and shape you need a resolution. If you don’t have specific resolution you’d like to warp to, you can use (as an approximation) the shape of the raster if you were to simply crop it without reprojecting:

src = rasterio.open(url)
dst_bounds_in_src_crs = rasterio.warp.transform_bounds(
    src.crs,
    dst_crs,
    left=src.bounds.left,
    right=src.bounds.right,
    top=src.bounds.top,
    bottom=src.bounds.bottom,
)
window_w_px = (dst_bounds_in_src_crs.right - dst_bounds_in_src_crs.left) / src.res[0]
window_h_px = (dst_bounds_in_src_crs.top - dst_bounds_in_src_crs.bottom) / src.res[1]
resolution = (bounds.right - bounds.left) / window_w_px, (bounds.top - bounds.bottom) / window_h_px

With this resolution, we can figure out the desired shape of the output. Rounding this output to an integer will modify the resolution slightly, but preserve the bounds.

width_approx = (bounds.right - bounds.left) / resolution[0]
height_approx = (bounds.top - bounds.bottom) / resolution[1]
dst_width = round(width_approx)
dst_height = round(height_approx)
dst_transform = rasterio.transform.from_bounds(
    west=bounds.left,
    south=bounds.bottom,
    east=bounds.right,
    north=bounds.top,
    width=dst_width,
    height=dst_height,
)

Now that we have the affine transform and shape, the warp can be done:

a_warped = np.zeros((dst_height, dst_width), dtype=np.float32)
a_warped[:] = np.nan
rasterio.warp.reproject(
    source=src.read(),
    destination=a_warped,
    src_crs=src.crs,
    dst_crs=dst_crs,
    src_transform=src.transform,
    dst_transform=dst_transform,
    dst_resolution=resolution,
)

Now you can write the new array back to disk:

profile = src.profile.copy()
profile.update(**kwargs)
profile["width"] = a_warped.shape[1]
profile["height"] = a_warped.shape[0]
profile["transform"] = dst_transform
profile["crs"] = dst_crs
with rasterio.open("warped.tif", "w", **profile) as f:
    f_write.write(a_warped, 1)