gis shp与shp相交实现

92 阅读 0 评论 61 点赞

我是靠谱客的博主魁梧小海豚，这篇文章主要介绍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
```
1
2
TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near
```
报错:TopologyException: Input geom 1 is invalid
尝试此种方法无果。
多边形自相交处理-selfIntersection
最后解决方案参考文章：
修复不合法的polygon
编码问题
最后综合的属性中，中文打印为乱码，因此将encoding设置为utf-8 即可解决

fiona使用过程中

复制代码

1
2
3
4
5
6
7
8
9
10
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()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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()