python gdal两个矢量shp叠加

在Python中使用GDAL库进行矢量叠加可以通过以下步骤实现：

1. 导入所需的库和模块：

import ogr
import osr
from osgeo import gdal

2. 打开两个矢量文件：

shp1 = ogr.Open("path/to/shp1.shp")
lyr1 = shp1.GetLayer()

shp2 = ogr.Open("path/to/shp2.shp")
lyr2 = shp2.GetLayer()

3. 创建一个新的矢量数据集和图层来存储叠加结果：

driver = ogr.GetDriverByName("ESRI Shapefile")
output = driver.CreateDataSource("path/to/output.shp")
output_lyr = output.CreateLayer("output", geom_type=ogr.wkbPolygon)

4. 复制shp1的属性表结构到输出图层中：

lyr1_def = lyr1.GetLayerDefn()
for i in range(lyr1_def.GetFieldCount()):
    output_lyr.CreateField(lyr1_def.GetFieldDefn(i))

5. 逐个读取shp1的要素，并使用Intersect方法与shp2进行叠加：

for feat1 in lyr1:
    geom1 = feat1.GetGeometryRef()

    # 筛选与shp1要素相交的shp2要素
    lyr2.SetSpatialFilter(geom1)

    for feat2 in lyr2:
        geom2 = feat2.GetGeometryRef()
        
        # 判断两个要素是否相交
        if geom1.Intersect(geom2):
            # 创建新的要素
            output_feat = ogr.Feature(output_lyr.GetLayerDefn())
            
            # 复制shp1的属性值到输出要素中
            for i in range(lyr1_def.GetFieldCount()):
                output_feat.SetField(i, feat1.GetField(i))

            # 将shp1和shp2的几何形状进行叠加
            output_feat.SetGeometry(geom1.Intersection(geom2))
            
            # 将新要素添加到输出图层
            output_lyr.CreateFeature(output_feat)
            
            # 删除shp2筛选器
            lyr2.SetSpatialFilter(None)
            break

6. 关闭打开的矢量文件和输出数据集：

shp1 = None
shp2 = None
output = None

这样，你就可以通过使用GDAL库的API在Python中实现两个矢量shp的叠加操作了。