遥感矢量文件坐标系转换到WGS1984

ytkz2025-02-192025-02-19

今天分享一个坐标系转换到wgs1984的脚本，并且分析此脚本生成的结果、和arcgis、QGIS的结果的差别，以此验证该脚本的可靠性。首先那么我说下为什么写这个脚本。

1

现有一个遥感矢量文件，它的坐标系是投影坐标，如下图。

需要把它的坐标系转为地理坐标wgs1984。

我知道arcgis或者QGis可以实现这个功能，但是这样的操作，是做不到批量处理的，所以有必要自己写代码实现这个功能，从而加深对矢量文件处理的理解。

2

下面的代码，是一个使用GDAL/OGR库的Python脚本，用来将输入的Shapefile文件转换到目标坐标系，默认是WGS84（EPSG:4326）

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2025/2/19 20:05
# @File : 矢量转为wgs1984坐标系.py
'''

'''
# !/usr/bin/env python
from osgeo import ogr, osr
import os


def reproject_shapefile(input_shp, output_shp, target_epsg=4326):
    """
    将输入的 shp 文件转换为目标坐标系（默认为WGS84, EPSG:4326），并保存到新的 shp 文件中。

    参数：
      input_shp   : str, 输入的 shp 文件路径（非WGS84坐标系）
      output_shp  : str, 输出的 shp 文件路径（转换后为WGS84坐标系）
      target_epsg : int, 目标坐标系EPSG编号，默认为4326（WGS84）
    """
    # 获取ESRI Shapefile驱动
    driver = ogr.GetDriverByName("ESRI Shapefile")

    # 打开输入的 shp 文件（只读模式）
    in_ds = driver.Open(input_shp, 0)
    if in_ds is None:
        raise FileNotFoundError(f"无法打开输入文件：{input_shp}")

    in_layer = in_ds.GetLayer()
    in_srs = in_layer.GetSpatialRef()  # 获取输入文件的空间参考

    # 设置输入空间参考的坐标轴映射策略为传统GIS顺序（经度，纬度）
    in_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    # 定义目标空间参考（WGS84）
    target_srs = osr.SpatialReference()
    target_srs.ImportFromEPSG(target_epsg)
    # 设置目标空间参考的坐标轴映射策略为传统GIS顺序
    target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    # 构建坐标转换对象：从输入坐标系转换到目标坐标系
    transform = osr.CoordinateTransformation(in_srs, target_srs)

    # 如果输出文件已存在，则删除它
    if os.path.exists(output_shp):
        driver.DeleteDataSource(output_shp)

    # 创建输出数据集
    out_ds = driver.CreateDataSource(output_shp)
    if out_ds is None:
        raise Exception(f"无法创建输出文件：{output_shp}")

    # 创建新的图层，采用与输入图层相同的几何类型，空间参考设为目标坐标系
    out_layer = out_ds.CreateLayer("reprojected", srs=target_srs, geom_type=in_layer.GetGeomType())

    # 将输入图层的字段结构复制到输出图层
    in_layer_defn = in_layer.GetLayerDefn()
    for i in range(in_layer_defn.GetFieldCount()):
        field_defn = in_layer_defn.GetFieldDefn(i)
        out_layer.CreateField(field_defn)
    out_layer_defn = out_layer.GetLayerDefn()

    # 遍历输入图层的所有要素，对其几何进行坐标转换后写入输出图层
    for in_feature in in_layer:
        geom = in_feature.GetGeometryRef()
        # 对几何进行转换
        geom.Transform(transform)

        # 创建新的要素，并设置几何和属性
        out_feature = ogr.Feature(out_layer_defn)
        out_feature.SetGeometry(geom)
        for i in range(out_layer_defn.GetFieldCount()):
            field_name = out_layer_defn.GetFieldDefn(i).GetNameRef()
            out_feature.SetField(field_name, in_feature.GetField(field_name))

        out_layer.CreateFeature(out_feature)
        out_feature = None  # 释放资源

    # 清理关闭数据集
    in_ds = None
    out_ds = None
    print(f"转换完成，输出文件保存在：{output_shp}")



if __name__ == "__main__":
    # 设置输入和输出 shp 文件路径
    input_shp = r'F:\vector.shp' # 替换为实际的输入文件路径
    output_shp = r'F:\vector1_wgs84.shp'  # 替换为你希望保存的输出文件路径
    reproject_shapefile(input_shp, output_shp)

3

这段代码的本质是执行矢量数据（Shapefile）的坐标重投影，将几何对象的坐标值从原始坐标系转换到目标坐标系（如WGS84）。

它不仅仅是修改元数据中的EPSG编码，而是对几何坐标进行数学变换，确保数据在新坐标系下位置准确。

代码关键步骤：

坐标转换对象：transform = osr.CoordinateTransformation(in_srs, target_srs) 定义了数学转换规则。
几何变换：geom.Transform(transform) 对每个几何对象应用转换，修改其坐标值。
输出文件：新图层使用目标坐标系元数据，并存储转换后的坐标。

其实，以上的代码还能支持对矢量文件的不同坐标系的转换，只需要修改target_epsg参数。

目前target_epsg参数是4326，即WGS1984坐标系。这个需要你提前查询你要转换的坐标系的epsg数值是什么。