【PythonGIS】矢量数据投影转换(坐标转换)

news2024/11/28 22:42:12

之前跟大家分享过面矢量数据投影转换和点矢量数据投影转换，但博主在日常工作的过程中发现之前分享的面矢量数据投影转换有时候会出现错误，或者转换后的效果不好。再一次偶然的过程中发现了新的坐标转换（投影转换）函数，所以今天和大家分享一下，大家在进行坐标转换时可以两种方法都试试，哪个效果好用哪个。

一、导入库

from osgeo import gdal, ogr, osr

二、坐标转换函数

这里最重要的就是gdal.VectorTranslate这个函数，这个函数包含了很多参数，代码中已给出相关注释，大家如果还有什么问题，可以直接去查看下官方文档。

此外我这里的目标坐标系是使用EPSG编码定义的，之前分享的文章中展示过从已有文件中提取坐标系，大家可以去看下然后放入代码即可。这里给大家分享两个查询EPSG编码的网站。

mygeodata.cloudhttps://mygeodata.cloud/cs2cs/

epsg.oihttps://epsg.io/

def Projection_Transform(shape_path, save_path, epsg):
    ogr.RegisterAll()
    # 注册所有驱动
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 定义编码类型为UTF8，防止字段乱码
    ds = ogr.Open(shape_path)
    layer = ds.GetLayer()
    # 打开需要转换的矢量数据，获取图层
    src_proj = layer.GetSpatialRef()
    # 获取其源坐标信息
    layer_name = layer.GetName()
    # 获取图层名称
    ds.Destroy()
    dst_proj = osr.SpatialReference()
    dst_proj.ImportFromEPSG(epsg)
    # 定义目标投影的空间信息
    options = gdal.VectorTranslateOptions(
        format="ESRI Shapefile",
        # 输出格式
        accessMode=None,
        srcSRS=src_proj,
        # 源坐标系
        dstSRS=dst_proj,
        # 目标坐标系
        reproject=True,
        selectFields=None,
        # 保留全部字段
        layerName=layer_name,
        # 图层名称
        geometryType="POLYGON",
        # 数据类型
        dim="XY"
    )
    gdal.VectorTranslate(save_path, srcDS=shape_path, options=options)

三、完整代码

# -*- coding: utf-8 -*-
"""
@Time ： 2023/8/22 17:59
@Auth ： RS迷途小书童
@File ：Vector Data Projection Transform.py
@IDE ：PyCharm
@Purpose ：矢量数据投影转换/坐标转换
"""
from osgeo import gdal, ogr, osr


def Projection_Transform(shape_path, save_path, epsg):
    ogr.RegisterAll()
    # 注册所有驱动
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 定义编码类型为UTF8，防止字段乱码
    ds = ogr.Open(shape_path)
    layer = ds.GetLayer()
    # 打开需要转换的矢量数据，获取图层
    src_proj = layer.GetSpatialRef()
    # 获取其源坐标信息
    layer_name = layer.GetName()
    # 获取图层名称
    ds.Destroy()
    dst_proj = osr.SpatialReference()
    dst_proj.ImportFromEPSG(epsg)
    # 定义目标投影的空间信息
    options = gdal.VectorTranslateOptions(
        format="ESRI Shapefile",
        # 输出格式
        accessMode=None,
        srcSRS=src_proj,
        # 源坐标系
        dstSRS=dst_proj,
        # 目标坐标系
        reproject=True,
        selectFields=None,
        # 保留全部字段
        layerName=layer_name,
        # 图层名称
        geometryType="POLYGON",
        # 数据类型
        dim="XY"
    )
    gdal.VectorTranslate(save_path, srcDS=shape_path, options=options)


if __name__ == "__main__":
    in_path = "B:/Personal/九段西.shp"
    out_path = "B:/Personal/九段西1.shp"
    in_path.replace("/", "\\")
    Projection_Transform(in_path, out_path, epsg=4550)
    # 输入文件，输出文件，目标坐标系