矢量数据投影转换-CFANZ编程社区

矢量数据投影转换

作者：阿振

案例说明

接着上一篇博文中，我们得到了WGS84坐标系下的中国省区图，而我们一般中国地图中使用的是割圆锥投影。

由于我国位于中纬度地区，中国地图和分省地图经常采用割圆锥投影，中国地图的中央经线常位于东经105度，两条标准纬线分别为北纬25度和北纬47度，而各省的参数可根据地理位置和轮廓形状初步加以判定。

在SpatialReference中查到我们一般使用的中国地图投影为：http://spatialreference.org/ref/sr-org/8657

PROJ4格式的定义为：+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

使用该投影，我们祖国雄鸡才会变得雄赳赳气昂昂，更好地展现我们神州大地的风采。

方法介绍

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

使用命令工具及其对应的命令行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()
使用GDAL提供的基本API进行实现
如果要自己利用基本API函数实现的话，基本思路如下：

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

代码实现

调用gdal.VectorTranslate()命令行工具的包装函数实现：

from osgeo import gdal
import os
os.environ['SHAPE_ENCODING'] = "utf-8"


src_file = 'China.shp'
dst_file = 'China_Reprojected.shp'

# 使用命令行API转换
# 输出数据投影定义，参考资料：http://spatialreference.org/ref/sr-org/8657
srs_def = """+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 """
gdal.VectorTranslate(dst_file, src_file, dstSRS=srs_def, reproject=True)

src_ds = ogr.Open(src_file)
src_layer = src_ds.GetLayer(0)
src_srs = src_layer.GetSpatialRef()  # 输入数据投影

调用基本API函数实现

from osgeo import ogr
from osgeo import osr
import os
os.environ['SHAPE_ENCODING'] = "utf-8"


src_file = 'China.shp'
dst_file = 'China_Reprojected.shp'

src_ds = ogr.Open(src_file)
src_layer = src_ds.GetLayer(0)
src_srs = src_layer.GetSpatialRef()  # 输入数据投影

# 输出数据投影定义，参考资料：http://spatialreference.org/ref/sr-org/8657
srs_def = """+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 """
dst_srs = osr.SpatialReference()
dst_srs.ImportFromProj4(srs_def)

# 创建转换对象
ctx = osr.CoordinateTransformation(src_srs, dst_srs)

# 创建输出文件
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(dst_file)
dst_layer = dst_ds.CreateLayer('province', dst_srs, ogr.wkbPolygon)

# 给输出文件图层添加属性定义
layer_def = src_layer.GetLayerDefn()
for i in range(layer_def.GetFieldCount()):
    field_def = layer_def.GetFieldDefn(i)
    dst_layer.CreateField(field_def)

# 循环遍历源Shapefile中的几何体添加到目标文件中
src_feature = src_layer.GetNextFeature()
while src_feature:
    geometry = src_feature.GetGeometryRef()
    geometry.Transform(ctx)
    dst_feature = ogr.Feature(layer_def)
    dst_feature.SetGeometry(geometry)  # 设置Geometry
    # 依次设置属性值
    for i in range(layer_def.GetFieldCount()):
        field_def = layer_def.GetFieldDefn(i)
        field_name = field_def.GetName()
        dst_feature.SetField(field_name, src_feature.GetField(field_name))
    dst_layer.CreateFeature(dst_feature)
    dst_feature = None
    src_feature = None
    src_feature = src_layer.GetNextFeature()
dst_ds.FlushCache()

del src_ds
del dst_ds

# 创建Shapefile的prj投影文件
dst_srs.MorphToESRI()
(dst_path, dst_name) = os.path.split(dst_file)
with open(dst_path + os.pathsep + dst_name + '.prj', 'w') as