矢量数据投影转换

- 日理万妓 2022-05-23 06:22 254阅读 0赞

# 矢量数据投影转换 #

作者：阿振

邮箱：tanzhenyugis@163.com

博客：[https://blog.csdn.net/theonegis/article/details/80089375][https_blog.csdn.net_theonegis_article_details_80089375]

修改时间：2018-06-03

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

--------------------

## 案例说明 ##

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

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

在[SpatialReference][]中查到我们一般使用的中国地图投影为：[http://spatialreference.org/ref/sr-org/8657][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()`命令行工具的包装函数实现：

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函数实现

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())

[https_blog.csdn.net_theonegis_article_details_80089375]: https://blog.csdn.net/theonegis/article/details/80089375
[SpatialReference]: http://spatialreference.org
[http_spatialreference.org_ref_sr-org_8657]: http://spatialreference.org/ref/sr-org/8657