GDAL重投影：投影变换和投影坐标转地理坐标

做 Python GIS 自动化时，GDAL重投影经常出现在三个具体任务里：把一幅地方投影的 GeoTIFF 转成 WGS84，经由脚本把 Shapefile 或 GeoPackage 统一到项目坐标系，以及把工程坐标、米制坐标转换成经纬度用于 WebGIS 定位。很多错误不是 GDAL 命令写错，而是源数据投影没有查清、把“设置投影”误当成“投影变换”，或者经纬度轴顺序被反了。本文按可复用流程讲清楚投影信息检查、栅格重投影、矢量转换和单点坐标转换。

GDAL重投影先做三件事：查投影、确认坐标、再变换

不要拿到数据就直接运行转换命令。正确顺序是先检查源数据是否真的带有坐标参考系统，也就是 CRS；再确认坐标值的单位和范围是否符合这个 CRS；最后才执行投影转换。如果第一步判断错，后面的结果看起来可能能打开，但会整体偏移几百米、几公里，甚至跑到另一个国家。

例如一个栅格的左上角坐标是 395000, 3456000，这明显不是经纬度，更像投影坐标。如果它被错误标成 EPSG:4326，地图软件可能会把 395000 度当作经度处理，结果自然无法正常显示。此时要先找出真实的源 CRS，再决定是否执行GDAL重投影。

问题背景：为什么投影转换后数据会偏移

在实际项目里，投影偏移通常来自四类情况。第一类是源数据没有 CRS，例如某些 CAD 转出的 Shapefile 只有坐标值，没有 .prj 文件；第二类是 CRS 写错，例如明明是 CGCS2000 高斯投影，却被标成 WGS84；第三类是把米制投影坐标当成经纬度使用；第四类是 Python 里没有处理好 EPSG:4326 的轴顺序，导致经度和纬度位置互换。

因此，投影处理不是单纯“转一下坐标系”。对栅格来说，投影变换还涉及输出范围、像元大小、重采样算法和 NoData；对矢量来说，投影变换会改写每个几何点的坐标；对一个单独的 X、Y 坐标来说，只需要坐标转换，不应该把整幅数据重新采样一遍。

这也是很多初学者踩坑的原因：GDAL重投影、投影坐标转经纬度、写入投影信息，三件事看起来相似，但实际影响完全不同。

核心原理：CRS、坐标值和像元网格分别代表什么

CRS 是坐标参考系统，说明坐标值如何对应到地球表面。常见的 EPSG:4326 是地理坐标系，单位通常是度；UTM、高斯克吕格、Web Mercator 和地方平面坐标多数是投影坐标系，单位通常是米。投影坐标转地理坐标，本质上是把同一个地面位置从米制平面坐标表达成经纬度。

对栅格数据，GDAL 还需要地理变换参数，也就是 GeoTransform。它描述像素列号、行号如何换算成地图坐标。只改变 CRS 字段不会移动像元，也不会重新计算每个像元的位置；真正的栅格投影变换会创建一个新的网格，并把源像元重采样到目标 CRS 中。

对矢量数据，几何点、线、面的每个坐标点都会被转换。它不涉及像元重采样，但会涉及几何精度、字段保留、图层驱动和目标格式限制。理解这一区别后，再看 gdalwarp、ogr2ogr、gdaltransform 和 Python osr API，就不会混用工具。

标准操作步骤：从检查到输出复查

无论处理栅格、矢量还是坐标表，都可以按下面流程走。它能把投影问题拆成可检查的几个环节。

1. 读取源数据空间参考。先用 gdalinfo、ogrinfo 或 Python API 查看 CRS、范围、单位和元数据。
2. 判断任务类型。如果只是补缺失 CRS，用赋值工具；如果要改变整份数据的坐标系，用重投影工具；如果只处理一个点或一张坐标表，用坐标转换工具。
3. 确认源 CRS。缺 CRS 时不要猜，优先查数据生产说明、同区域控制点、项目坐标系统或已知底图。
4. 选择目标 CRS。WebGIS 展示常见目标是 EPSG:4326 或 EPSG:3857，工程分析则应按项目要求选择地方投影或统一平面坐标系。
5. 执行转换并保存新文件。不要直接覆盖原始数据，保留可回退的源文件和处理日志。
6. 复查输出结果。再次查看 CRS、范围、像元大小或几何范围，并叠加到底图或控制点上验证。

GDAL获取投影信息：先确认源数据的 CRS

GDAL获取投影信息是所有转换前的第一步。命令行最直接的方法是 gdalinfo。它会输出驱动、尺寸、坐标系、仿射变换、范围、波段和 NoData 等信息。只看文件名或图层名称不可靠，必须以数据内部记录和项目元数据为准。

gdalinfo source.tif

如果输出里能看到 Coordinate System is、ID["EPSG",xxxx]、角点坐标和像元大小，说明这份栅格至少包含可读取的空间参考。若坐标系为空、角点坐标像普通像素坐标，或者范围明显不合理，就不要直接转换，应先回到数据来源确认真实 CRS。

在 Python 中，可以读取投影 WKT、EPSG 代码和 GeoTransform。下面示例适合批量检查影像文件，尤其适合在数据入库或批处理前做质检。

from osgeo import gdal, osr

path = "source.tif"
ds = gdal.Open(path)
if ds is None:
    raise RuntimeError("无法打开数据")

projection = ds.GetProjection()
geotransform = ds.GetGeoTransform()

if not projection:
    print("没有读取到投影信息")
else:
    srs = osr.SpatialReference()
    srs.ImportFromWkt(projection)
    print("CRS 名称:", srs.GetName())
    print("Authority:", srs.GetAuthorityName(None), srs.GetAuthorityCode(None))
    print("GeoTransform:", geotransform)

如果 GetAuthorityCode(None) 输出为空，不一定代表数据没有 CRS。有些数据只有 WKT，没有明确 EPSG 编号；有些地方坐标系或工程坐标系本来就没有公开 EPSG。此时应保存 WKT，并结合数据提供方说明、控制点或已知底图进行核对。读取投影信息的目的不是只拿到一个 EPSG，而是确认源坐标的真实含义。

GDAL投影变换：栅格重投影用 gdalwarp

GDAL投影变换用于把数据从一个 CRS 转到另一个 CRS。栅格最常用的是 gdalwarp，因为它会重新计算输出网格、范围和像元值。下面示例把源栅格转为 EPSG:4326。如果源数据内部已经有正确 CRS，通常只需要指定目标 CRS。

gdalwarp -t_srs EPSG:4326 -r bilinear -dstnodata -9999 source.tif source_4326.tif

如果源文件缺少 CRS，但你已经从数据说明中确认它的真实 CRS，可以在转换时用 -s_srs 指定源 CRS。注意，这不是凭空猜测投影，而是在命令中告诉 GDAL 源坐标应该如何解释。

gdalwarp -s_srs EPSG:4547 -t_srs EPSG:4326 -r bilinear -dstnodata -9999 source.tif source_4326.tif

重采样算法要按数据类型选择。连续型栅格，例如 DEM、温度、降雨、遥感反射率，可以用 bilinear 或 cubic；分类栅格，例如土地利用、行政区编码、分类结果，应该用 near，否则类别值会被插值成不存在的新值。

Python 自动化中可以使用 gdal.Warp。它适合把多个文件接入批处理脚本，统一输出坐标系、压缩方式和 NoData。

from osgeo import gdal

gdal.Warp(
    "source_4326.tif",
    "source.tif",
    dstSRS="EPSG:4326",
    resampleAlg="bilinear",
    dstNodata=-9999,
    creationOptions=["TILED=YES", "COMPRESS=DEFLATE"]
)

做GDAL重投影后，一定要再次运行 gdalinfo 检查输出文件。重点看目标 CRS 是否正确、角点经纬度是否落在合理区域、像元大小是否符合预期、NoData 是否保留。能打开文件不代表结果正确。

矢量数据的投影变换用 ogr2ogr

矢量数据属于 OGR 处理范围，但在日常使用中也归入 GDAL 工具链。把 Shapefile、GeoPackage、GeoJSON 或 PostGIS 图层转到目标 CRS，常用 ogr2ogr。下面示例把地块数据转换到 EPSG:4326，并输出为 GeoPackage。

ogr2ogr -t_srs EPSG:4326 parcels_4326.gpkg parcels.gpkg

如果源矢量缺少坐标系，而你确认它实际是 EPSG:4547，可以同时指定源 CRS 和目标 CRS。

ogr2ogr -s_srs EPSG:4547 -t_srs EPSG:4326 parcels_4326.gpkg parcels_without_prj.shp

这里同样要区分两件事：-s_srs 是说明源坐标如何解释，-t_srs 是目标输出坐标系。只有指定目标 CRS，几何坐标才会被转换。如果只是给一个没有 .prj 的 Shapefile 补投影文件，不能当作真正的坐标变换。

GDAL投影坐标转换成地理坐标：单点和批量坐标怎么转

GDAL投影坐标转换成地理坐标适合处理单个点、坐标表、鼠标拾取坐标、栅格像元中心点等场景。它只转换坐标值，不重采样栅格，也不改写整份矢量数据。如果你只是想把 x=395000, y=3456000 转成经纬度，使用 gdaltransform 或 Python osr.CoordinateTransformation 更直接。

echo 395000 3456000 | gdaltransform -s_srs EPSG:4547 -t_srs EPSG:4326

在 Python 中，建议显式设置传统 GIS 轴顺序。这样输入投影坐标时按 x, y，输出地理坐标时按 longitude, latitude 理解，减少 EPSG:4326 轴顺序带来的混淆。

from osgeo import osr

src = osr.SpatialReference()
src.ImportFromEPSG(4547)

dst = osr.SpatialReference()
dst.ImportFromEPSG(4326)

src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

transform = osr.CoordinateTransformation(src, dst)

x, y = 395000, 3456000
lon, lat, z = transform.TransformPoint(x, y)
print(lon, lat)

如果要把栅格的像素行列号转换成经纬度，需要先用 GeoTransform 把像素坐标转成地图坐标，再把地图坐标转成地理坐标。下面示例计算某个像元中心点的经纬度。

from osgeo import gdal, osr

ds = gdal.Open("source.tif")
gt = ds.GetGeoTransform()

src = osr.SpatialReference()
src.ImportFromWkt(ds.GetProjection())
src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

dst = osr.SpatialReference()
dst.ImportFromEPSG(4326)
dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

transform = osr.CoordinateTransformation(src, dst)

col, row = 1200, 800
x = gt[0] + (col + 0.5) * gt[1] + (row + 0.5) * gt[2]
y = gt[3] + (col + 0.5) * gt[4] + (row + 0.5) * gt[5]

lon, lat, z = transform.TransformPoint(x, y)
print(lon, lat)

这类脚本很适合把影像上的检测结果、样本点、像元中心点或切片索引定位到 WebGIS 地图上。要注意的是，这种坐标转经纬度只解决坐标表达问题，不会修复源数据本身的投影错误。

设置投影和重投影不是一回事

这是使用 GDAL 时最重要的边界。设置投影是给数据写入或覆盖 CRS 元数据，坐标值不变；重投影是根据源 CRS 和目标 CRS 计算新的坐标值或新的栅格网格。把这两件事混用，是投影偏移的高频原因。

例如一幅 GeoTIFF 没有投影，但你确认它本来就是 EPSG:4547，可以用 gdal_translate -a_srs 写入投影信息。这个命令不会把米制坐标变成经纬度。

gdal_translate -a_srs EPSG:4547 source_no_srs.tif source_with_srs.tif

如果要把它真正转成 EPSG:4326，则需要继续使用 gdalwarp。前者是补标签，后者是计算新位置。判断自己的任务属于哪一种，是避免错误的关键。

常见坑点与排查方法

• 把错误 CRS 写进源数据。如果真实源 CRS 不确定，不要为了让软件“能显示”随便写 EPSG:4326。先找数据说明、同区域控制点或权威底图核对。
• 把 gdal_translate -a_srs 当作重投影。它只设置空间参考，不重新计算坐标。需要改变坐标系时，应使用 gdalwarp 或 ogr2ogr -t_srs。
• 经纬度顺序反了。Python 坐标转换中要特别注意轴顺序。面向常规 GIS 工作流时，可设置 OAMS_TRADITIONAL_GIS_ORDER，按经度、纬度处理输出。
• 分类栅格用了双线性重采样。土地利用、分类编码、掩膜栅格应使用 nearest neighbour，也就是 -r near。
• NoData 没有传递到输出。重投影后如果背景值参与插值，边界会出现异常颜色或无效像元扩散。应明确设置源和目标 NoData。
• 输出分辨率没有控制。目标 CRS 单位变化后，默认像元大小可能不符合项目需要。必要时用 -tr 或 -ts 控制输出分辨率或尺寸。
• 只看图层是否能打开。投影错误的数据也可能在软件里显示。应检查角点坐标、范围、叠加底图位置和已知点误差。

如果GDAL重投影后结果偏移，优先回查源 CRS 和轴顺序，而不是反复试不同目标 EPSG。投影问题通常不是“多试几次”能解决的，必须先把坐标语义确认清楚。

工具和方法对比

不同任务应选择不同工具。下面表格可以作为日常判断依据。

实践检查清单

每次处理投影转换前，可以按下面清单快速自检。它比直接套命令更可靠。

• 源数据是否能通过 gdalinfo 或脚本读取 CRS。
• 坐标值范围是否符合 CRS 单位，例如经纬度不应出现几十万这样的坐标值。
• 是否已经确认源 CRS，而不是凭显示效果猜测。
• 任务是补投影信息、整份数据重投影，还是单点坐标转换。
• 栅格数据是否选择了适合的重采样算法。
• NoData、输出分辨率、压缩和内部瓦片是否需要明确设置。
• Python 坐标转换是否处理了轴顺序。
• 输出结果是否再次检查 CRS、角点坐标、范围和底图叠加位置。

把这份清单固定到项目脚本里，能让GDAL重投影从一次性手工操作变成可追溯的数据处理流程。

FAQ：GDAL重投影常见问题

GDAL获取投影信息只看到 WKT，看不到 EPSG 怎么办？

这很常见。WKT 本身已经能描述 CRS，不一定必须有 EPSG 编号。可以先检查 WKT 里的投影名称、椭球、中央经线、单位和坐标范围，再结合数据来源说明确认。如果项目必须统一到 EPSG 编号，可以用 GIS 软件或权威坐标系统资料做匹配，但不要仅凭名称相似就替换。

GDAL投影变换和设置投影有什么区别？

这个过程会把坐标从源 CRS 计算到目标 CRS，栅格会生成新网格，矢量会改写几何坐标。设置投影只是给数据写入 CRS 元数据，坐标值不变。缺 CRS 但坐标本来正确时可以设置投影；需要从米制坐标转经纬度时，必须执行真正的转换。

GDAL投影坐标转换成地理坐标后为什么经纬度反了？

多半是轴顺序问题。EPSG:4326 在不同上下文中可能按纬度、经度或经度、纬度理解。面向常规 GIS 和 WebGIS 使用时，Python 脚本里建议设置 OAMS_TRADITIONAL_GIS_ORDER，并明确把输出解释为 longitude, latitude。

投影坐标转经纬度需要重投影整幅栅格吗？

不需要。如果只是要把某个 X、Y 坐标、像元中心点或样本点转成经纬度，用 gdaltransform 或 osr.CoordinateTransformation 就够了。只有当你需要输出一份新的栅格数据，并让整幅图处于目标 CRS 下，才需要用 gdalwarp。

重投影后图像边界变斜或范围变大正常吗？

正常。投影变换会把源数据边界投到新的坐标系中，矩形边界在目标 CRS 下可能不再保持原来的方向。GDAL 会为输出栅格计算能覆盖结果的外接范围，所以范围可能变大，边缘也可能出现 NoData 区域。

总结

GDAL重投影的关键不是记住某个命令，而是先判断任务类型：读取投影信息用于确认源数据坐标语义，整份数据投影转换用于把栅格或矢量转到目标 CRS，单点坐标转经纬度用于坐标表或像元中心点定位。只要把“查清源 CRS、区分设置投影和重投影、处理重采样与轴顺序、输出后复查范围”这四步做好，投影偏移问题会少很多，Python GIS 批处理也更容易复用到真实项目中。