How to extract a region of interest (ROI) from a raster with rasterio. This will maintain the original projection of the source data. Data are converted to Cloud Optimezed GeoTIFF (COG) format via rio_cogeo

rio-roi.py

import rasterio
from rasterio.mask import mask
import rio_cogeo
from shapely.geometry import box


def region(src, output_path, geom, geom_crs):
    region_geom = rasterio.warp.transform_geom(geom_crs, src.crs, geom)

    # Read the raster data within the region
    out_image, out_transform = mask(src, [region_geom], crop=True)

    # Update the metadata
    profile = src.meta.copy()
    profile.update({
        'height': out_image.shape[1],
        'width': out_image.shape[2],
        'transform': out_transform,
        'driver': 'GTiff',
    })
    profile.update({
        'blockxsize': 256, 
        'blockysize': 256,
        'compress': 'lzw', 
        'quality': 90, # this may not be a valid rio option
        'tiled': True,  # Issue: this doesn't build overviews... use rio-cogeo
        'bigtiff': 'IF_SAFER',
    })
    print(profile)

    with rasterio.open(output_path, 'w', **profile) as dst:
        dst.write(out_image)
    # Use rio-cogeo to produce overviews properly
    rio_cogeo.cog_translate(output_path, output_path, profile, use_cog_driver=True)


# Define the input raster file path
input_raster_path = 'path/to/source.tif'

# Define the output raster file path
output_raster_path = 'output_raster.tif'

# Define the coordinates of the region in EPSG:4326 CRS
# Specify the bounding box coordinates as (min_longitude, min_latitude, max_longitude, max_latitude)
region_crs = 'EPSG:4326'
# region_coordinates = (-74.1, 40.6, -73.9, 40.8)
region_coordinates = (-109.1382, 36.9148, -102.0081, 41.0462)
region_geom = box(*region_coordinates)