gis shp与shp相交实现

场景

工作要解决的问题：

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

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

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

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

碰到的问题

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

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

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

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

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

4. 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 进行传递成功
   

5. 有序字典的合并
   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()

最后

以上就是魁梧小海豚最近收集整理的关于gis shp与shp相交实现的全部内容，更多相关gis内容请搜索靠谱客的其他文章。