【Arcpy】快速计算一个要素类中两个要素相交面积百分比

小鱼儿 2023-07-24 13:17 11阅读 0赞

## \# 问题所在 ##

一个要素类（也就是一个图层内）中有两类，比如`地块`和`建筑`。

我们想要统计一个数据，即`地块`内的`建筑面积总数`，这可能存在两种情况：

*  **情形1**：建筑完全在地块内，面积可直接加到建筑面积总数上
 *  **情形2**：建筑和地块边缘相交，我们要判断建筑和建筑相交面积站建筑自身面积的百分比，若相交面积大于建筑面积50%，则加到建筑面积上，否则忽略

通过分析可以发现，判定两个要素是否相交并且计算他们相交面积占建筑要素的百分比是难点所在。

按正常思路，我们需要做双重循环，逐个遍历地块和建筑要素，判断他们是包含关系还是相交关系，在此基础上再计算相交面积百分比。问题是，如果建筑要素数量极大，比如超过10000个，这种双重循环将耗费大量时间。

另一个问题是，如何在一个要素类中求得两个要素之间的相交部分？

## \# 解决方案 ##

通过多次试验探索，总结出一条可行的思路：

1.  基于空间位置选择要素：使用`arcpy.SelectLayerByLocation_management`，以地块要素图形为界限，选择其周围的建筑要素；
2.  判定要素空间关系：步骤1已经选出要素类中所有和地块存在重叠、包含、相交等关系的要素，可以使用`geometry`的`contains`、`equals`方法做条件判定，过滤出相交情形；
3.  相交要素插入临时要素类，做相交分析：把相交的的两个要素分别插入（`insertRow`）事先建立的临时要素类，使用相交分析(`Intersect_analysis`)，得到相交要素类，获得相交要素类和建筑要素类的面积并计算比例；相交分析完成后，分别清空临时要素类并删去相交要素类；
4.  面积汇总及更新地块要素。

## \# 参考代码 ##

注意：运行代码前在`工作空间.gdb`中提前准备好`地块建筑_合并`要素类，在`临时数据.gdb`提前准备好`big`和`small`要素类：

*  **地块建筑\_合并**：地块和建筑两个图层合并后的图层，需要额外添加`Category`和`Sum_Area`字段。`Category`用来区分地块和建筑，其中，`Category=0`的是地块要素，`Category=1`的为建筑要素；`Sum_Area`用来记录该地块内所有建筑面积之和。
 *  **big**：一个空图层，用来临时保存所计算的地块要素的图层，用于计算和建筑要素类的交集，可以直接复制`地块建筑_合并`要素类，然后删去其中的所有要素。
 *  **small**：一个空图层，用来临时保存所计算的建筑要素的图层，用于计算和地块要素类的交集，可以直接复制`地块建筑_合并`要素类，然后删去其中的所有要素。

# -*- coding: utf-8 -*-
    import arcpy
    
    fc_path = unicode("./工作空间.gdb/地块建筑_合并", "utf-8")
    big_path = unicode("./临时数据.gdb/big", "utf-8")
    small_path = unicode("./临时数据.gdb/small", "utf-8")
    intersect_path = unicode("./临时数据.gdb/intersect", "utf-8")
    
    arcpy.DeleteRows_management(big_path)
    arcpy.DeleteRows_management(small_path)
    arcpy.Delete_management(intersect_path)
    
    bigs = arcpy.da.UpdateCursor(fc_path,("SHAPE@", "Sum_Area"), "Category=0 And Sum_Area is Null")
    big_insert = arcpy.InsertCursor(big_path)
    small_insert = arcpy.InsertCursor(small_path)
    
    # 需要更新的地块要素，Category=0的为地块要素
    bigs = arcpy.da.UpdateCursor(fc_path,("SHAPE@", "Sum_Area"), "Category=0")
    
    # 使用原数据创建图层，用于空间选择操作
    arcpy.MakeFeatureLayer_management(fc_path, "lyr")
    
    num = 0
    
    for big in bigs:
        print "########## next shape ##########"
        print num
        big_geometry = big[0]
        sum_area = 0
    	
    	# 基于地块空间位置选择其等同、包含、相交的所有建筑要素
        arcpy.SelectLayerByLocation_management("lyr",'INTERSECT', big_geometry)
        smalls = arcpy.da.SearchCursor('lyr', ("SHAPE@", "Shape_Area"), "Category=1")
    
        for small in smalls:
            small_geometry = small[0]
            small_area = small[1]
    
            if big_geometry.equals(small_geometry):
                continue
    
            elif big_geometry.contains(small_geometry):
                sum_area += small_area 
            
            else:
            	# 将相交的一个地块要素和一个建筑要素分别插入临时要素类，用于计算相交要素类
                big_row = big_insert.newRow()
                big_row.shape = big_geometry
                big_insert.insertRow(big_row)
    
                small_row = big_insert.newRow()
                small_row.shape = small_geometry
                small_insert.insertRow(small_row)
                
                # try:
                # arcpy.Intersect_analysis([big_path, small_path], intersect_path)
                # except:
                # print arcpy.GetMessages().encode('gbk')
    
                arcpy.Intersect_analysis([big_path, small_path], intersect_path)
    
    			# 如果没有相交内容（由于select和intersect的相交标准不一致导致相交分析没有结果） 
                try:
                    intersect_area = arcpy.da.SearchCursor(intersect_path, ("Shape_Area")).next()[0]
                except:
                    arcpy.DeleteRows_management(big_path)
                    arcpy.DeleteRows_management(small_path)
                    arcpy.Delete_management(intersect_path)
                    continue
                 
                intersect_percentage = intersect_area / small_area
                print "%s%" % str(round(intersect_percentage * 100, 2))
                if intersect_percentage >= 0.5:
                    sum_area += small_area
    
                # 清空临时要素类，删除相交要素类
                arcpy.DeleteRows_management(big_path)
                arcpy.DeleteRows_management(small_path)
                arcpy.Delete_management(intersect_path)
    
    
        print "sum area: ", sum_area
        big[1] = sum_area
        bigs.updateRow(big)
        num += 1