Python+GDAL实现面矢量数据叠加

需要使用python和GDAL实现矢量数据的相交叠加计算，空间上的叠加分析，本来是一个非常简单的功能。但是在网上却只能找到c++版的相关代码。Python的代码几乎找不到，找一个发现是把两个shape文件中的图形逐个取出来相交，效率十分低下。最后根据C++的代码，实现了Python版本的调用GDAL实现两个矢量shape文件叠加分析。

环境

• python3.9
• GDAL 3.4.3

代码

from osgeo import ogr, osr
from osgeo import gdal
import os


def calculated_proportion(in_path1, in_path2, out_path):
    """
    面图层相交，属性字段保留，两个shp文件必须是同一个坐标系
    :param in_path1: 主shp文件
    :param in_path2: 用来相交的shp文件
    :param out_path: 输出结果路径
    :return:
    """
    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = driver.Open(in_path1, 1)
    layer_1 = dataSource.GetLayer()
    crs = layer_1.GetSpatialRef()
    epsg1 = crs.GetAttrValue('AUTHORITY', 1)

    RdataSource = driver.Open(in_path2, 1)
    layer_2 = RdataSource.GetLayer()
    crs = layer_2.GetSpatialRef()
    epsg2 = crs.GetAttrValue('AUTHORITY', 1)
    if epsg1 != epsg2:
        print("坐标系不一致")
    else:
        print("坐标系一致")

    # 新建DataSource，Layer
    out_ds = driver.CreateDataSource(out_path)
    out_lyr = out_ds.CreateLayer("out", layer_1.GetSpatialRef(), layer_1.GetGeomType())
    # 遍历原始的Shapefile文件给每个Geometry做Buffer操作
    # current_union = layer[0].Clone()
    print('the length of layer:', len(layer_1))
    if len(layer_1) == 0:
        return

    p = ["SKIP_FAILURES=YES", "PROMOTE_TO_MULTI=YES", "INPUT_PREFIX=1", "METHOD_PREFIX=2"]
    # gdal.TermProgress_nocb 控制台输出进度
    layer_1.Intersection(layer_2, out_lyr, p, gdal.TermProgress_nocb)
    out_lyr.SyncToDisk()
    del dataSource, RdataSource, out_ds

if __name__ == '__main__':
    shp_path = "E:/测试数据/DK4107262017.shp"
    shp_path2 = "E:/测试数据/zzjgProj.shp"
    out_path = "E:/测试数据/out.shp"
    calculated_proportion(shp_path, shp_path2, out_path)

效果

注意

网上有人提到，相交后，结果是空的。我也遇到了这个问题，最后发现是两个叠加图层的坐标系不一致就会出现这个问题。