July 11, 2018

矢量数据投影转换

矢量数据投影转换

作者:阿振

邮箱:tanzhenyugis@163.com

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

修改时间:2018-06-03

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


案例说明

接着上一篇博文中,我们得到了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库,我们有两种方法进行矢量数据的重投影:

  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()命令行工具的包装函数实现:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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函数实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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())