Rasterio坐标转换?Reproject怎么写?

wowwwai GIS研习社 · 工具流程与项目排障

很多同学在处理栅格数据时会卡在“Rasterio坐标转换?Reproject怎么写?”这个问题上：明明知道输入影像是一个坐标系，输出需要变成另一个坐标系，但不知道 CRS、transform、宽高和重采样参数该怎么一起配置。本文用一个可复用的 Python 示例，把 Rasterio 坐标转换、Rasterio reproject 写法和常见报错排查串起来。

引言：Rasterio坐标转换到底在转换什么

在 Rasterio 里，栅格坐标转换通常不是简单改一个 EPSG 编号，而是要把影像像元重新放到新的空间参考系统中。这个过程一般叫栅格重投影，在 Rasterio 中主要通过 rasterio.warp.reproject 完成。

初学者最容易混淆三件事：

• CRS：坐标参考系统，例如 EPSG:4326、EPSG:3857、EPSG:4490。
• transform：仿射变换参数，用来描述像元行列号和真实地图坐标之间的关系。
• reproject：把源栅格按目标 CRS 和目标 transform 重新采样到新网格。

所以，Rasterio reproject 写法的重点不是只填 dst_crs，还要正确计算 dst_transform、dst_width、dst_height，否则输出影像可能位置偏移、范围错误，或者像元大小不符合预期。

背景：什么时候需要用 Rasterio reproject

以下场景都属于典型的 Rasterio坐标转换需求：

• 把 WGS84 经纬度影像从 EPSG:4326 转为 WebGIS 常用的 EPSG:3857。
• 把遥感影像转为项目使用的投影坐标系，例如 UTM、CGCS2000 高斯投影。
• 让多期栅格数据统一到同一个 CRS，方便栅格叠加分析。
• 把 DEM、土地利用、NDVI 等栅格统一坐标系后再做裁剪、统计或制图。
• 修复“矢量和栅格叠不上”的问题，但前提是原始 CRS 本身是正确的。

需要注意：如果影像本来没有 CRS，或者 CRS 标错了，直接 reproject 并不能自动修正。Rasterio 只能按照你提供的源 CRS 去计算坐标转换，源信息错了，结果也会跟着错。

原理：Rasterio坐标转换的四个关键参数

理解 Rasterio reproject 之前，先看一个栅格数据在空间中的定位方式。一个 GeoTIFF 影像至少依赖这几类信息：

• 像元数组：真实的栅格值，例如高程、分类编码、反射率。
• 源 CRS：源影像所在坐标系，对应 src.crs。
• 源 transform：源影像的行列号到地图坐标的换算关系，对应 src.transform。
• 宽高：源影像的列数和行数，对应 src.width、src.height。

当你把栅格从一个坐标系转到另一个坐标系时，Rasterio 需要回答三个问题：

1. 目标坐标系是什么？也就是 dst_crs。
2. 目标影像覆盖多大范围？通常由源影像范围转换得到。
3. 目标影像像元如何排列？也就是 dst_transform、dst_width、dst_height。

Rasterio 提供的 calculate_default_transform 就是用来根据源 CRS、目标 CRS、源范围和源宽高，计算目标 transform、目标宽度和目标高度。

from rasterio.warp import calculate_default_transform

dst_transform, dst_width, dst_height = calculate_default_transform(
    src_crs,
    dst_crs,
    src_width,
    src_height,
    *src_bounds
)

有了这些参数后，再调用 reproject 把每个波段从源网格重采样到目标网格。

步骤：Rasterio reproject 怎么写完整代码

步骤一：安装 Rasterio

建议在独立的 Python 环境中安装 Rasterio。使用 conda 通常更省心，因为 GDAL 相关依赖会一起处理。

conda install -c conda-forge rasterio

如果使用 pip，也可以安装：

pip install rasterio

如果你在 Windows 上遇到 GDAL 依赖问题，优先考虑 conda-forge 环境。

步骤二：查看输入影像 CRS 和 transform

在做 Rasterio坐标转换之前，先确认输入影像是否带有正确 CRS。

import rasterio

input_tif = "input.tif"

with rasterio.open(input_tif) as src:
    print("CRS:", src.crs)
    print("Transform:", src.transform)
    print("Width:", src.width)
    print("Height:", src.height)
    print("Bounds:", src.bounds)
    print("Count:", src.count)

正常情况下，src.crs 应该能输出类似下面的结果：

CRS: EPSG:4326

如果输出是 None，说明影像没有写入坐标系信息。此时不要急着 reproject，要先确认原始数据到底是什么坐标系。

步骤三：把 GeoTIFF 从 EPSG:4326 转为 EPSG:3857

下面是一段完整、可复用的 Rasterio reproject 代码。它会读取输入 GeoTIFF，计算目标网格，并输出一个新的投影坐标系 GeoTIFF。

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

input_tif = "input.tif"
output_tif = "output_epsg3857.tif"
dst_crs = "EPSG:3857"

with rasterio.open(input_tif) as src:
    transform, width, height = calculate_default_transform(
        src.crs,
        dst_crs,
        src.width,
        src.height,
        *src.bounds
    )

    kwargs = src.meta.copy()
    kwargs.update({
        "crs": dst_crs,
        "transform": transform,
        "width": width,
        "height": height
    })

    with rasterio.open(output_tif, "w", **kwargs) as dst:
        for band_index in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, band_index),
                destination=rasterio.band(dst, band_index),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest
            )

这段代码中最关键的是：

• calculate_default_transform：计算目标影像的 transform、width、height。
• kwargs.update：更新输出 GeoTIFF 的空间参考和尺寸信息。
• reproject：逐波段执行重投影。
• Resampling.nearest：使用最近邻重采样，适合分类栅格。

步骤四：根据数据类型选择重采样方法

Rasterio reproject 必须指定重采样方式。不同数据类型应选择不同方法。

数据类型	推荐重采样	原因
土地利用、行政区编码、分类结果	Resampling.nearest	保持类别编码不被插值改变
DEM、高程、温度、降水、NDVI	Resampling.bilinear	连续变量更平滑，适合插值
需要更平滑的连续栅格	Resampling.cubic	视觉效果更平滑，但计算更慢
汇总到更粗分辨率	Resampling.average	适合连续变量降尺度平均

例如，对 DEM 做坐标转换时，可以这样写：

resampling=Resampling.bilinear

而对土地利用分类数据，通常不要使用双线性或三次卷积，否则原本的类别值可能被插值成无意义的小数。

步骤五：验证输出结果是否正确

完成 Rasterio坐标转换后，不要只看文件是否生成。建议至少检查四项：

import rasterio

output_tif = "output_epsg3857.tif"

with rasterio.open(output_tif) as dst:
    print("CRS:", dst.crs)
    print("Transform:", dst.transform)
    print("Width:", dst.width)
    print("Height:", dst.height)
    print("Bounds:", dst.bounds)
    print("Count:", dst.count)

重点看：

• CRS 是否变成目标坐标系。
• Bounds 是否处于目标坐标系的合理范围。
• Width 和 Height 是否没有异常变得极大或极小。
• 把输出影像加载到 QGIS 或 ArcGIS Pro 中，是否能和同 CRS 的底图、矢量数据正确叠加。

常见坑：Rasterio坐标转换失败或结果错位的原因

1. 只修改 CRS，没有真正重投影

有些人会直接修改元数据：

kwargs.update({"crs": "EPSG:3857"})

如果只改 CRS，而不调用 reproject，这相当于“给影像换标签”，像元位置并没有重新计算。结果通常是图层严重错位。

判断原则：如果源坐标系和目标坐标系不同，通常需要 reproject；如果只是原文件缺少 CRS，但坐标值本来就是正确的，则需要补写 CRS，而不是重投影。

2. 输入影像 src.crs 是 None

如果 src.crs 是 None，下面这类代码会失败：

calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)

解决思路是先确认原始影像真实坐标系。如果你确认它是 EPSG:4326，可以在读取时显式指定源 CRS 参与 reproject：

src_crs = "EPSG:4326"

但这必须建立在你确定数据来源和坐标定义的基础上，不能靠猜。

3. 重采样方法选错

分类栅格使用 bilinear 或 cubic 后，可能出现 1.3、2.7 这类不存在的类别值。连续栅格使用 nearest 则可能出现明显锯齿。Rasterio reproject 写法没有绝对统一模板，重采样方法必须根据数据含义选择。

4. nodata 没有处理好

如果输入影像有 NoData 值，建议在输出元数据中保留，并在 reproject 中明确设置：

reproject(
    source=rasterio.band(src, band_index),
    destination=rasterio.band(dst, band_index),
    src_transform=src.transform,
    src_crs=src.crs,
    dst_transform=transform,
    dst_crs=dst_crs,
    src_nodata=src.nodata,
    dst_nodata=src.nodata,
    resampling=Resampling.nearest
)

否则边界区域、无效区域可能被错误填充，影响后续统计。

5. 输出文件过大

从经纬度坐标系转为投影坐标系时，如果目标分辨率计算异常，输出宽高可能非常大。遇到这种情况，先打印 width、height 和 transform，确认目标像元大小是否合理。

print(width, height)
print(transform)

如果需要控制输出分辨率，可以在更高级的流程中手动设置目标 transform 和输出尺寸，但这需要你明确目标像元大小和目标范围。

方法比较：Rasterio reproject、GDAL Warp 和 QGIS 重投影

方法	适合场景	优点	注意点
Rasterio reproject	Python 自动化处理、多文件批处理、和 NumPy 分析结合	代码可控，适合构建脚本和数据处理流水线	需要理解 CRS、transform、重采样参数
GDAL Warp	命令行批处理、服务端数据转换	成熟稳定，参数丰富，性能好	命令参数较多，新手容易漏选项
QGIS 栅格重投影	少量数据、交互式检查、教学演示	界面直观，便于查看结果	自动化能力不如 Python 脚本
ArcGIS Pro Project Raster	ArcGIS 工作流、企业项目制图与分析	与地理处理工具链集成好	需要注意环境参数和重采样设置

如果你只是偶尔转换一个影像，QGIS 或 ArcGIS Pro 更直观。如果你要处理几十个或几百个栅格文件，Rasterio坐标转换脚本更适合复用和自动化。

检查清单：写 Rasterio reproject 前后逐项确认

写 Rasterio reproject 代码时，可以按下面的清单检查，能避免大多数坐标转换问题。

• 输入影像是否能正常打开？
• src.crs 是否存在且正确？
• src.transform 是否合理？
• src.bounds 是否符合源坐标系的数值范围？
• 目标 CRS 是否写成正确的 EPSG 编码或 CRS 字符串？
• 是否使用 calculate_default_transform 计算目标 transform、width、height？
• 输出元数据是否更新了 crs、transform、width、height？
• 是否逐波段调用 reproject？
• 重采样方法是否符合数据类型？
• NoData 值是否正确传递？
• 输出结果是否在 QGIS 或 ArcGIS Pro 中叠加检查过？

FAQ：Rasterio坐标转换常见问题

Q1：Rasterio reproject 和修改 crs 有什么区别？

修改 CRS 只是改变数据的坐标系标签，不改变像元位置。Rasterio reproject 会根据源 CRS 和目标 CRS 重新计算像元位置，并把数据重采样到目标网格。如果两个坐标系不同，通常应该使用 reproject。

Q2：为什么 Rasterio坐标转换后影像位置偏了？

常见原因有三个：源 CRS 本身标错、只改了 CRS 没有重投影、目标 transform 或宽高设置错误。建议先打印 src.crs、src.bounds、dst.crs、dst.bounds，再用 QGIS 加载源数据和结果数据进行对比。

Q3：Rasterio reproject 可以处理多波段影像吗？

可以。常见写法是遍历 range(1, src.count + 1)，对每个波段调用一次 reproject。本文示例已经包含多波段处理逻辑。

Q4：分类栅格坐标转换应该用什么 resampling？

分类栅格建议使用 Resampling.nearest。例如土地利用类型、分类编码、行政区栅格都不应该被插值成小数。连续栅格如 DEM、温度、降水可以考虑 bilinear 或 cubic。

Q5：输入影像没有 CRS，能不能直接指定 EPSG 后 reproject？

只有在你明确知道原始影像真实坐标系时才可以。比如数据说明文档明确写着是 WGS84 经纬度坐标，那么可以把源 CRS 设为 EPSG:4326 后再转换。不要凭肉眼猜 CRS，否则输出结果很可能错误。

Q6：Rasterio坐标转换后像元大小为什么变了？

从一个 CRS 转到另一个 CRS 时，像元形状和单位会发生变化。例如从经纬度的“度”转到 Web Mercator 的“米”，Rasterio 会根据源范围和目标坐标系计算默认目标分辨率。因此输出像元大小变化是正常现象，但需要检查是否符合分析需求。

结论：Rasterio reproject 的正确写法要同时处理 CRS、transform 和重采样

解决“Rasterio坐标转换?Reproject怎么写?”这个问题，关键不是背一段固定代码，而是理解完整流程：读取源影像空间信息，确定目标 CRS，用 calculate_default_transform 计算目标网格，再用 reproject 按正确重采样方式逐波段写出。

实际项目中，建议把本文代码作为基础模板，再根据数据类型补充 NoData、压缩、分辨率控制和批处理逻辑。只要确认源 CRS 正确、目标参数合理、重采样方法匹配数据含义，Rasterio坐标转换就会稳定很多。