扫二维码与项目经理沟通
我们在微信上24小时期待你的声音
解答本文疑问/技术咨询/运营咨询/技术建议/互联网交流
需要使用python和GDAL实现矢量数据的相交叠加计算,空间上的叠加分析,本来是一个非常简单的功能。但是在网上却只能找到c++版的相关代码。Python的代码几乎找不到,找一个发现是把两个shape文件中的图形逐个取出来相交,效率十分低下。最后根据C++的代码,实现了Python版本的调用GDAL实现两个矢量shape文件叠加分析。
成都创新互联公司是专业的墨玉网站建设公司,墨玉接单;提供网站设计制作、成都网站设计,网页设计,网站设计,建网站,PHP网站建设等专业做网站服务;采用PHP框架,可快速的进行墨玉网站开发网页制作和功能扩展;专业做搜索引擎喜爱的网站,专业的做网站团队,希望更多企业前来合作!环境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元起,快前往官网查看详情吧
我们在微信上24小时期待你的声音
解答本文疑问/技术咨询/运营咨询/技术建议/互联网交流