# 矢量数据投影转换

作者：阿振

邮箱：<tanzhenyugis@163.com>

博客：<https://blog.csdn.net/theonegis/article/details/80089375>

修改时间：2018-06-03

声明：本文为博主原创文章，转载请注明原文出处

## 案例说明

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

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

在[SpatialReference](http://spatialreference.org)中查到我们一般使用的中国地图投影为：<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库，我们有两种方法进行矢量数据的重投影：

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及其属性写入目标文件

## 代码实现

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

```python
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()  # 输入数据投影
```

1. 调用基本API函数实现

```python
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 f:
    f.write(dst_srs.ExportToWkt())
```


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://theonegis.gitbook.io/geopy/gdal-kong-jian-shu-ju-chu-li/gdal-shu-ju-ji-ben-cao-zuo/shi-liang-shu-ju-tou-ying-zhuan-huan.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.