gis shp与shp相交实现

场景

工作要解决的问题：

项目根据卫星影像和算法得到同一区域两个时间的变化情况，可以得到一个包含变化图斑的.shp文件，但是该区域中不同区域的性质不同，需要进行区别处理。

依据是甲方提供的一个区域分割的图斑.shp文件。

最后，需要输出的是甲方每个图斑中真实的变化情况(即前shp文件和后shp文件求交集)，还要将两个shp中的图斑属性进行综合。

最先参考文章 OGR Layer Intersection 之后(最终使用方案) Intersect Shapefiles using Shapely 之后还参考 Fiona用户手册

碰到的问题

首先碰到的是两个shp无交集，但是软件做有，最后发现是两个shp的坐标参考系不同

在进行图斑相交处理过程中会碰到

TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near

报错:TopologyException: Input geom 1 is invalid 尝试此种方法无果。 多边形自相交处理-selfIntersection 最后解决方案参考文章： 修复不合法的polygon

编码问题 最后综合的属性中，中文打印为乱码，因此将encoding设置为utf-8 即可解决

fiona使用过程中

with fiona.open( 'oriented-ccw.shp', 'w', crs=source.crs, # 该参数报错 driver=source.driver, schema=sink_schema, ) as sink: # 值后尝试此种 sink = fiona.open('/tmp/foo.shp', 'w', **source.meta) 最后传参使用 crs_wkt 进行传递成功

有序字典的合并 Merge 两个 OrderedDict 到一个 OrderedDict 中的可用和不可用方法

代码

import time import fiona import rtree from shapely.geometry import shape, mapping def intersect_shp_shp(manager_shp_path,input_shp_path,output_shp_path,min_area=None): """ manager_shp_path : 第一shp文件路径 该项目中是 林地小班 input_shp_path : 第二shp文件路径 该项目中是 算法产生的变化图斑 output_shp_path : 输出shp文件路径 min_area : 用于图斑面积筛选(面积小于该值的图斑，将会pass，也提高算法执行时间) return None """ start_=time.time() with fiona.open(manager_shp_path, 'r',encoding='utf-8') as layer1: with fiona.open(input_shp_path, 'r',encoding='utf-8') as layer2: if layer1.crs['init']!=layer2.crs['init']: raise Exception('shapefile的坐标系不同，无法进行计算！请先转换坐标系！') # 合并两shp的schema schema3 = layer2.schema.copy() schema3['properties'].update(layer1.schema['properties']) # 新增面积属性 schema3['properties']['area'] = 'float' # schema3['geometry']='Polygon' with fiona.open(output_shp_path, mode='w', schema=schema3,driver='ESRI Shapefile',crs_wkt=layer2.meta['crs_wkt'],encoding='utf-8') as layer3: # 建立rtree索引 print('开始建立索引......') index = rtree.index.Index() manager_num=0 for feat1 in layer1: fid = int(feat1['id']) geom1 = shape(feat1['geometry']) index.insert(fid, geom1.bounds) manager_num+=1 print('林地小班数量共有：%d 个，建立索引共耗时%.2f s'%(manager_num,time.time()-start_)) # 执行相交运算 intersect_execute(layer1,layer2,layer3,index,min_area) print('完成本次运算共耗时%.2f s'%(time.time()-start_,)) def intersect_execute(layer1,layer2,layer3,index,min_area): print('开始进行交集运算......') result_num,polygon_num,small_num=0,0,0 for feat2 in layer2: polygon_num+=1 geom2 = shape(feat2['geometry']) if min_area and geom2.area<min_area: small_num+=1 continue # 检测合法性，并进行合法化(主要针对 有洞的polygon) if not geom2.is_valid: geom2=geom2.buffer(0) for fid in list(index.intersection(geom2.bounds)): feat1 = layer1[fid] geom1 = shape(feat1['geometry']) if not geom1.is_valid: geom1=geom1.buffer(0) if geom1.intersects(geom2): # 合并属性 props = feat2['properties'].copy() props.update(feat1['properties']) intersect=geom1.intersection(geom2) if min_area and intersect.area<min_area: continue # 给area属性赋值 props['area']=intersect.area try: layer3.write({ 'properties': props, 'geometry': mapping(intersect) }) result_num+=1 except Exception as e: print(e) print('输入图斑数量 %d 个，小面积的图斑有 %d 个，得到符合要求的图斑共： %d 个'%(polygon_num,small_num,result_num,)) if __name__ == "__main__": intersect_shp_shp()