如何用Rasterio从非EPSG:4326坐标系数据中裁剪指定范围？

我想裁剪/加载大型数据集的一部分，但裁剪范围采用EPSG:4326坐标系，而数据集本身并非该坐标系。尝试多种计算方式后，仍无法获取正确的裁剪窗口。

当前代码 #

import rasterio as rio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds

def getWindow(dst_crs, dst_transform):
    south, west, north, east = S2_box.split(", ")
    newbounds = transform_bounds(
        rio.CRS.from_epsg(4326),
        dst_crs,
        float(west),
        float(south),
        float(east),
        float(north),
    )
    window = from_bounds(*newbounds, transform=dst_transform)
    return window

def S2_TCI(ds, name):
    """Creates Sentinel 2 true color image (TCI)"""
    name = f"{name}-TCI"
    print(name)
    sds = rio.open(ds.GetSubDatasets()[c.DS_TCI][0])
    profile = sds.profile
    bands = sds.read(
        [c.BAND_RED, c.BAND_GRN, c.BAND_BLU],
        window=getWindow(sds.crs, sds.transform),
    )
    writeTiffRGB(bands, profile, name)

问题排查与修正 #

核心问题点 #

• 坐标顺序传递误差：transform_bounds要求输入参数顺序为(src_crs, dst_crs, left, bottom, right, top)，你拆分S2_box得到的south, west, north, east对应bottom, left, top, right，虽已调整顺序，但需确保完全匹配函数要求的left(west), bottom(south), right(east), top(north)逻辑。
• 未处理窗口越界：转换后的边界可能超出数据集实际范围，导致读取空数据或报错，需将生成的窗口与数据集有效范围取交集。

修正后的代码 #

import rasterio as rio
from rasterio.windows import from_bounds, intersection
from rasterio.warp import transform_bounds

def getWindow(dst_crs, dst_transform, dst_width, dst_height):
    # 拆分EPSG:4326格式的边界：south, west, north, east
    south, west, north, east = S2_box.split(", ")
    # 转换坐标系，严格遵循left(west), bottom(south), right(east), top(north)的顺序
    new_bounds = transform_bounds(
        rio.CRS.from_epsg(4326),
        dst_crs,
        float(west),
        float(south),
        float(east),
        float(north),
    )
    # 基于转换后的边界生成裁剪窗口
    window = from_bounds(*new_bounds, transform=dst_transform)
    # 确保窗口不超出数据集的有效范围
    valid_window = intersection(window, rio.windows.Window(0, 0, dst_width, dst_height))
    return valid_window

def S2_TCI(ds, name):
    """生成哨兵2号真彩色图像（TCI）"""
    name = f"{name}-TCI"
    print(name)
    sds = rio.open(ds.GetSubDatasets()[c.DS_TCI][0])
    profile = sds.profile.copy()
    # 获取经过有效性校验的裁剪窗口
    window = getWindow(sds.crs, sds.transform, sds.width, sds.height)
    # 读取指定波段的裁剪后数据
    bands = sds.read(
        [c.BAND_RED, c.BAND_GRN, c.BAND_BLU],
        window=window
    )
    # 更新输出图像的元数据，匹配裁剪后的尺寸和坐标变换
    profile.update(
        width=window.width,
        height=window.height,
        transform=rio.windows.transform(window, sds.transform)
    )
    writeTiffRGB(bands, profile, name)

额外注意事项 #

• 确保S2_box的格式为south, west, north, east，例如"40.1, -105.2, 40.3, -105.0"，无多余空格或格式错误。
• 调用writeTiffRGB前必须更新profile的宽高和变换参数，否则输出图像会出现坐标偏移或尺寸不符的问题。

内容的提问来源于stack exchange，提问作者Stefan Gofferje