矢量坐标转换

image-20231222153133081

需要把非84坐标的矢量转换为WGS84坐标的矢量。

方法介绍

跟栅格数据投影转换一样，使用GDAL库，我们有两种方法进行矢量数据的重投影：

1.使用命令工具及其对应的命令行API接口进行转换（简单，准确，实践中一定要用这种方法）

GDAL提供了ogr2ogr命令行工具进行矢量数据投影转换，命令如下：ogr2ogr -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=30 +lon_0=105 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs " China_Projected.shp China.shp

-t_srs选项制定输出数据投影，当然可以是ESPG，也可以是PROJ4或者OGC WKT格式的投影定义都OK

GDAL对该命令的封装的C/C++函数是GDALVectorTranslate(),Python中是gdal.VectorTranslate()

2.使用GDAL提供的基本API进行实现。如果要自己利用基本API函数实现的话，基本思路如下：

  • 利用osgeo.ogr.Driver.CreateDataSource()创建输出数据
  • 根据源文件创建目标文件的属性字段定义
  • 利用osgeo.osr.CoordinateTransformation对象将源文件中的Geometry对象转为目标文件中的Geometry对象（其实质是进行不同投影系统下空间几何体的坐标转换）
  • 遍历源文件，依次将所有几何体的Geometry及其属性写入目标文件

代码

代码如下：

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2023/12/22 19:12 
# @File : transfer_project.py
from osgeo import ogr, gdal
from osgeo import osr
import os

def VectorTranslate(
        shapeFilePath,
        saveFolderPath,
        format="GeoJSON",
        accessMode=None,
        dstSrsESPG=4326,
        selectFields=None,
        geometryType="POLYGON",
        dim="XY",
):
    """
    转换矢量文件，包括坐标系，名称，格式，字段，类型，纬度等。
    :param shapeFilePath: 要转换的矢量文件
    :param saveFolderPath: 生成矢量文件保存目录
    :param format: 矢量文件格式，强烈建议不要使用ESRI Shapefile格式。
    :param accessMode:None代表creation,'update','append','overwrite'
    :param dstSrsESPG: 目标坐标系EPSG代码，4326是wgs84地理坐标系
    :param selectFields: 需要保留的字段列表如果都保留，则为None
    :param geometryType: 几何类型,"POLYGON","POINT"。。。
    :param dim: 新矢量文件坐标纬度,建议查阅官方API。
    :return:
    """
    ogr.RegisterAll()
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
    data = ogr.Open(shapeFilePath)
    layer = data.GetLayer()
    spatial = layer.GetSpatialRef()
    layerName = layer.GetName()
    data.Destroy()
    dstSRS = osr.SpatialReference()
    dstSRS.ImportFromEPSG(int(dstSrsESPG))
    if format == "GeoJSON":
        destDataName = layerName + ".geojson"
        destDataPath = os.path.join(saveFolderPath, destDataName)
    elif format == "ESRI Shapefile":
        destDataName = os.path.join(saveFolderPath, layerName)
        flag = os.path.exists(destDataName)
        os.makedirs(destDataName) if not flag else None
        destDataPath = os.path.join(destDataName, layerName + ".shp")
    else:
        print("不支持该格式！")
        return
    options = gdal.VectorTranslateOptions(
        format=format,
        accessMode=accessMode,
        srcSRS=spatial,
        dstSRS=dstSRS,
        reproject=True,
        selectFields=selectFields,
        layerName=layerName,
        geometryType=geometryType,
        dim=dim
    )
    gdal.VectorTranslate(
        destDataPath,
        srcDS=shapeFilePath,
        options=options
    )
    return destDataPath

if __name__ == '__main__':
    shapeFilePath = r'D:\china1820.shp'
    saveFolderPath = r'D:\china'
    VectorTranslate(
        shapeFilePath,
        saveFolderPath,
        format="ESRI Shapefile",
        accessMode=None,
        dstSrsESPG=4326,
        selectFields=None,
        geometryType="POLYGON",
        dim="XY",
    )

image-20231222153111974