ArcGis 计算Boyce-Clark形状指数（python）

早八起得来

已于 2024-04-08 14:09:38 修改

阅读量979

点赞数 10

分类专栏： arcgis 文章标签： arcgis python 空间计算

于 2023-12-09 20:40:14 首次发布

本文链接：https://blog.csdn.net/yuan_hou_/article/details/134900009

版权

arcgis 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

虾神原创，公众号\知乎：虾神说D

引用大神文章：新版白话空间统计（37）：几何形状规则度与Boyce-Clark半径形状指数 - 知乎

许多编程新手直接看大神的代码文章看不懂，这里对大神的文章进行梳理手把手教你如何计算Boyce-Clark形状指数。

点击“python”打开编辑窗口

把数据拖入map中

1. 导入arcpy包

import arcpy

设置数据处理环境，包括工作空间和允许数据覆盖，定义空间参考

#第一行路径根据情况进行修改

#第三行坐标根据情况进行修改（WKID码）

arcpy.env.workspace = r"E:\Other\1\data\PythonDemo\PythonDemo.gdb"
arcpy.env.overwriteOutput = True
sr = arcpy.SpatialReference(4326)

创建结果图层，并且添加需要的字段

pline = arcpy.CreateFeatureclass_management(arcpy.env.workspace, "sline", "POLYLINE", spatial_reference = sr)
arcpy.AddField_management(pline, "name", "TEXT")
bcipnt = arcpy.CreateFeatureclass_management(arcpy.env.workspace, "bci", "POINT", spatial_reference = sr)
arcpy.AddField_management(bcipnt, "name", "TEXT")
arcpy.AddField_management(bcipnt, "bci", "FLOAT")

点击回车生成两个要素。

进行下一步前先把地名的字段名称改为“Name”

2. 定义各种方法

子方法1 ：获取从中心点到extent的最大范围的方法

def maxdist(pnt, ext):
    point = arcpy.Point(ext.XMin, ext.YMin)
    pg1 = arcpy.PointGeometry(point, sr)
    dist1 = pnt.angleAndDistanceTo(pg1)[1]
    point = arcpy.Point(ext.XMax, ext.YMax)
    pg2 = arcpy.PointGeometry(point, sr)
    dist2 = pnt.angleAndDistanceTo(pg2)[1]
    return max(dist1, dist2)

子方法2：根据角度和距离，创建一条射线

def createInPnt(cent, dist, angle):
    adpnt = cent.pointFromAngleAndDistance(angle, dist)
    line = arcpy.Polyline(arcpy.Array([cent.centroid, adpnt.centroid]), sr)
    return line

子方法3：如果有多个交点，有两种处理方法，一种是取短线，一种是取长线

#选一个方法不要两个都输入（我没有试过）区别可以看开头引用的文章

def minPoint(cent, pnt):
    p = pnt.centroid
    dmin = cent.angleAndDistanceTo(arcpy.PointGeometry(p, sr))[1]
    for i in range(pnt.pointCount):
        dist1 = cent.angleAndDistanceTo(arcpy.PointGeometry(pnt[i], sr))[1]
        if dmin > dist1:
            dmin = dist1
            p = pnt[i]
    return (dmin, p)

def maxPoint(cent, pnt):
    dmax = 0
    p = pnt.centroid
    for i in range(pnt.pointCount):
        dist1 = cent.angleAndDistanceTo(arcpy.PointGeometry(pnt[i], sr))[1]
        if dmax < dist1:
            dmax = dist1
            p = pnt[i]
    return (dmax, p)

子方法4：将生成的结果的线，写入到数据中

def createPline(name, line):
    cursor = arcpy.da.InsertCursor("sline", ["SHAPE@", "name"])
    for l in line:
        row = (l, name)
        cursor.insertRow(row)
    del cursor

子方法5：计算BC指数

def culabcIndex(distlist, n):
    distsum = sum(distlist)
    bci = 0.0
    for d in distlist:
        bci += math.fabs(((d / distsum) * 100) - (100.0 / n))
    return bci

子方法6: 主功能

def bcindex(pline, name, poly, n):
    cent = arcpy.PointGeometry(poly.centroid, sr)
    angle = 360.0 / n
    ext = poly.extent
    dist = maxdist(cent, ext)
    distlist = []
    linelist = []
    for i in range(n):
        a = i * angle
        line = createInPnt(cent, dist, a)
        pnt = line.intersect(poly, 1)
        mp = maxPoint(cent, pnt)
        linelist.append(arcpy.Polyline(arcpy.Array([cent.centroid, mp[1]]), sr))
        distlist.append(mp[0])
    bci = culabcIndex(distlist, n)
    createPline(name, linelist)
    return bci

main函数

#第三行"cn"中的cn根据情况改名称（你需要计算的要素的名称）

fields = ['SHAPE@', 'name']
biclist = []
geo = [row for row in arcpy.da.SearchCursor("cn", fields)]
for g in geo:
    bci = bcindex(pline, g[1], g[0], 36)
    biclist.append((g[0].centroid, g[1], bci))
bcicursor = arcpy.da.InsertCursor(bcipnt, ["SHAPE@", "name", "bci"])
for b in biclist:
    print(b)
    row = b
    bcicursor.insertRow(row)
del bcicursor

点击回车自动计算