GDAL 实现投影转换

前言

❝

在GIS开发中，对于数据源首先要确定的第一件事就是坐标系统，这就涉及到坐标转换处理问题。而经常遇到的便是Shapefile数据的投影转换，如何高效、准确的将源数据坐标系转换到目标坐标系是我们需要研究解决的问题。

在之前的文章中讲了如何使用GDAL或者ogr2ogr工具将txt以及csv文本数据转换为Shp格式，本篇教程在之前一系列文章的基础上讲解如何使用GDAL 实现投影转换。

• GDAL 简介
• GDAL 下载安装
• GDAL 开发起步
• GDAL 实现 GIS 数据读取转换（全）

如果你还没有看过，建议从以上内容开始。

1. 开发环境

本文使用如下开发环境，以供参考。

2. 数据准备

如下是本文选取的全国省级行政区Shp数据：

部分景点数据：

3. 投影形式

参考网站：[https://epsg.io/4490](https://epsg.io/4490)

坐标定义信息具有多种格式，可以选择符合通用标准的OGC格式，下面EPSG编码4490，即CGCS2000坐标系为例进行展示。

• OGC WKT

GEOGCS["China Geodetic Coordinate System 2000",
DATUM["China_2000",
SPHEROID["CGCS2000",6378137,298.257222101,
AUTHORITY["EPSG","1024"]],
AUTHORITY["EPSG","1043"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4490"]]

• ESRI WKT

GEOGCS["GCS_China_Geodetic_Coordinate_System_2000",
DATUM["D_China_2000",
SPHEROID["CGCS2000",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]

• PROJ.4

+proj=longlat +ellps=GRS80 +no_defs +type=crs

• Proj4js

proj4.defs("EPSG:4490","+proj=longlat +ellps=GRS80 +no_defs +type=crs");

• GeoServer

4490=GEOGCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4490"]]

4. 获取图层坐标系统

使用GetSpatialRef方法可以获取图层坐标参考，若图层缺少投影信息，则返回None。

# 获取图层及坐标系统
sourceLayer = shpDs.GetLayer(0)
geomType = sourceLayer.GetGeomType()
sourceSrs = sourceLayer.GetSpatialRef()

也可以使用Geometry对象方法GetSpatialReference获取。

# 获取坐标系统
geom = feature.GetGeometryRef()
sourceSrs = geom.GetSpatialReference()

5. 图层投影转换

5.1. 导入依赖

Shp作为一种矢量数据格式，可以使用矢量库OGR进行处理，用于打开数据源和获取图层。还需要引入osr模块用于坐标定义以及os模块用于判断文件数据路径。

from osgeo import ogr,osr
import os

5.2. 获取数据

定义一个方法LayerProject用于图层投影转换，该方法接收两个参数，一个源数据文件路径sourcePath，一个投影转换文件数据路径projectPath。

"""
说明：GDAL 投影转换
参数：
-sourcePath：源文件Shp数据路径
-projectPath：投影转换图层数据路径
"""
def LayerProject(sourcePath,projectPath):

按照老规矩添加数据驱动，使用checkFilePath方法检查文件路径是否存在，使用checkDriver判断数据驱动是否正常，之后获取图层及坐标系统。

# 检查文件是否存在
checkFilePath(sourcePath)
checkFilePath(projectPath)

# 获取数据驱动
shpDriver = ogr.GetDriverByName("ESRI Shapefile")

# 检查数据驱动是否正常
checkDriver(shpDriver)

# 获取数据源
shpDs = shpDriver.Open(sourcePath)

# 获取源图层及坐标信息
sourceLayer = shpDs.GetLayer(0)
geomType = sourceLayer.GetGeomType()
sourceSrs = sourceLayer.GetSpatialRef()

# 获取源数据结构
featureDefn = sourceLayer.GetLayerDefn()
fieldCount = featureDefn.GetFieldCount()

文件和数据驱动检查方法定义如下。

"""
说明：检查文件路径是否正常
参数：
-filePath：文件数据路径
"""
def checkFilePath(filePath):
ifos.path.exists(filePath):
print(f"{filePath} 文件数据路径存在")
else:
print(f"{filePath} 文件数据路径不存在，请检查！")

"""
说明：检查数据驱动是否正常
"""
def checkDriver(driver):
ifdriver is None:
print("数据驱动不可用")
returnFalse

5.3. 创建投影

使用osr.SpatialReference()创建空间参考，然后将投影信息导入到坐标系统。可以使用如下多种方法：

- ImportFromEPSG(SpatialReference self, int arg)
- ImportFromESRI(SpatialReference self, char ** ppszInput)
- ImportFromProj4(SpatialReference self, char * ppszInput)
- ImportFromUSGS(SpatialReference self, long proj_code, long zone=0, double [15] argin=0, long datum_code=0)
- ImportFromWkt(SpatialReference self, char ** ppszInput)

其中个人觉得最简单方便的还是ImportFromEPSG方法。

5.4. 创建投影图层

使用CreateDataSource方法创建投影数据源和图层，并根据源数据复制要素。创建图层方法CreateLayer第二个参数用于指定数据坐标系，坐标系统可以使用ImportFromEPSG方法定义，只需传入一个EPSG编码。

"""
投影转换操作
"""
# 创建投影数据源
prjDs = shpDriver.CreateDataSource(projectPath)

# 投影坐标对象
srs = osr.SpatialReference()
srs.ImportFromEPSG(4522)

# 创建投影图层
prjLayer = prjDs.CreateLayer("prj_layer",srs,geomType)

# 添加属性结构
foriinrange(fieldCount):
fieldDefn = featureDefn.GetFieldDefn(i)
prjLayer.CreateField(fieldDefn)

# 写入要素
forfeatureinsourceLayer:
prjLayer.CreateFeature(feature)

5.5. 导出投影

可以使用以下方法导出投影信息。

- ExportToPrettyWkt(SpatialReference self, int simplify=0)
- ExportToProj4(SpatialReference self)
- ExportToUSGS(SpatialReference self)
- ExportToWkt(SpatialReference self, char ** options=None)
- ExportToXML(SpatialReference self, char const * dialect="")

5.6. 几何投影

使用CoordinateTransformation方法定义转换信息，第一个参数为源数据坐标系，第二个参数为目标投影坐标系，然后调用Geometry对象方法Transform进行坐标转换。

# 源数据坐标参考
sourceSrs = osr.SpatialReference()
sourceSrs.ImportFromEPSG(4326)
print("源数据坐标系名称：",sourceSrs.GetName())

# 添加几何对象
geom = feature.GetGeometryRef()
print("之前的坐标：",geom.GetX(),geom.GetY())

# 几何投影
targetSrs = osr.SpatialReference()
targetSrs.ImportFromEPSG(4522)
print("目标数据坐标系名称：",targetSrs.GetName())

# 坐标转换
coordsTransform = osr.CoordinateTransformation(sourceSrs,targetSrs)
geom.Transform(coordsTransform)

print("之后的坐标：",geom.GetX(),geom.GetY())

坐标系输出信息显示如下。

5.7. 创建投影文件

对于缺少.prj文件的数据可以使用空间参考方法MorphToESRI()修改坐标信息，使用ExportToWkt()方法导出坐标参考后将其写入投影文件中。

创建一个方法CreatePrjFile用于创建投影文件。

"""
说明：创建.prj文件
参数：
-shpPath： Shp文件路径
"""
def CreatePrjFile(shpPath):
prjSrs = osr.SpatialReference()
prjSrs.ImportFromEPSG(4522)

prjSrs.MorphToESRI()

fileName = os.path.splitext(shpPath)[0]
prjFile = fileName +".prj"

with open(prjFile,"w") as f:
f.write(prjSrs.ExportToWkt())
print(f"成功创建投影文件: {prjFile}")

6. 注意事项

在windows开发环境中同时安装GDAL与PostGIS，其中投影库PROJ的环境变量指向PostGIS的安装路径，在运行GDAL程序时，涉及到要素、几何与投影操作时会导致异常。具体意思为GDAL不支持PostGIS插件中的投影库版本，需要更换投影库或者升级版本。

解决办法为修改PROJ的环境变量到GDAL支持的版本或者在GDAL程序开头添加以下代码：

os.environ['PROJ_LIB'] = r'D:\Programs\Python\Python311\Libsite-packages\osgeo\data\proj'