需要使用python和GDAL实现矢量数据的相交叠加计算，空间上的叠加分析，本来是一个非常简单的功能。但是在网上却只能找到c++版的相关代码。Python的代码几乎找不到，找一个发现是把两个shape文件中的图形逐个取出来相交，效率十分低下。最后根据C++的代码，实现了Python版本的调用GDAL实现两个矢量shape文件叠加分析。
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)
网上有人提到，相交后，结果是空的。我也遇到了这个问题，最后发现是两个叠加图层的坐标系不一致就会出现这个问题。
你是否还在寻找稳定的海外服务器提供商？创新互联www.cdcxhl.cn海外机房具备T级流量清洗系统配攻击溯源，准确流量调度确保服务器高可用性，企业级服务器适合批量采购，新人活动首月15元起，快前往官网查看详情吧
