Rasterio坐标转换?Reproject怎么写?

你写的Reproject为啥总报错？可能是坐标系理解错了

上周一位读者在后台留言：‘Dr. Gis，我用Rasterio做reproject，结果图像全黑，控制台还跳出CRS mismatch警告……’——这太典型了。我在参与全国耕地遥感监测项目时，也栽过这个坑：明明代码照着文档敲，输出却像被泼了墨。问题根源往往不在语法，而在你对‘坐标转换’本质的理解偏差。

坐标转换不是简单的像素搬家，而是空间关系的重构。就像把橘子皮强行摊平成世界地图——必然有拉伸、撕裂或重叠。

坐标系的本质：地球的“身份证”与“穿衣法则”

很多人以为CRS（Coordinate Reference System）只是个数字ID，其实它包含两层核心信息：

• 地理坐标系（Geographic CRS）：定义地球形状（椭球体）和原点（如WGS84），相当于地球的“身份证”。
• 投影坐标系（Projected CRS）：规定如何把曲面展开成平面（如UTM、Albers），是地球的“穿衣法则”——不同地区得穿不同尺码才合身。

当你调用rasterio.warp.reproject时，本质上是在执行“换装手术”：先脱掉旧投影（逆变换回地理坐标），再套上新投影（正变换到目标平面）。中间任何一步坐标系描述错误，都会导致像素“迷路”。

三步写出零报错的Reproject代码

以将Landsat影像从WGS84转为Albers等积圆锥投影为例，这是国土调查的黄金标准：

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

with rasterio.open('input.tif') as src:
    # 步骤1：计算目标栅格的仿射变换矩阵和尺寸
    dst_crs = 'EPSG:5070'  # Albers USGS
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    
    # 步骤2：创建目标文件并写入元数据
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })
    
    with rasterio.open('output.tif', 'w', **kwargs) as dst:
        # 步骤3：逐波段重投影（关键！避免内存爆炸）
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.bilinear  # 连续数据用bilinear，分类数据用nearest
            )

这段代码我优化过上百次——最易错的是calculate_default_transform的参数顺序，以及忘记更新kwargs中的宽高。曾经有实习生漏了这步，导致输出只有左上角1/4区域有数据，其余全是NoData。

避坑指南：三个高频雷区

终极心法：先验证，再生产

永远别相信第一次跑通的代码。我的习惯是：

1. 用QGIS叠加原图和重投影结果，肉眼检查边界是否对齐；
2. 抽样读取几个已知坐标的像素值（比如城市中心点），看数值是否合理；
3. 打印dst.bounds和dst.crs，确认元数据已更新。

记住：坐标转换的终极目标不是让代码不报错，而是让空间位置在地球上真实归位。你现在手头的项目，遇到过哪些奇葩的坐标转换问题？评论区留下你的血泪史，我来帮你诊断！