【ArcGIS】arcpy实现获取polygon的MBR角度

目标:想要获得建筑物的方向。通常以多边形的最小外包矩形的方向表示,即计算最小外包矩形及其角度。
思路:多边形→最小外包矩形→角度计算→赋值回多边形

7871333-c44347ec36bc8f27.png
最小外包矩形

方法与问题:ArcGIS中提供了计算最小外包矩形的工具;按照下图model builder思路进行可以实现,但ArcGIS未给出计算角度的工具(也可能我没找到,有知道的请告知,谢谢!),因此考虑使用arcpy实现此功能(为了省事儿直接全过程都arcpy了)

7871333-9a1120d292140a62.png
MBR Tool

7871333-e76e55ea69a81384.png
思路步骤

全部代码如下,作为个人记录,也分享给大家。。。

# -*- coding: utf-8 -*----
# Import arcpy module
import arcpy
import numpy as np

def distance(p1,p2):
    return ((p1.X-p2.X)**2+(p1.Y-p2.Y)**2)**0.5

def angle(p1,p2):
    dx=p1.X-p2.X
    dy=p1.Y-p2.Y
    return np.arctan(dy/dx)/np.pi*180

arcpy.env.overwriteOutput=True

# InputAddress
inputFeatureClass = arcpy.GetParameterAsText(0)
outputFeatureClass=arcpy.GetParameterAsText(1)

#输出文件
if arcpy.Exists(outputFeatureClass):
    arcpy.Delete_management(outputFeatureClass)
arcpy.CopyFeatures_management(inputFeatureClass,outputFeatureClass)
#外包矩形
if arcpy.Exists(arcpy.env.workspace+'/'+"gg"):
    arcpy.Delete_management(arcpy.env.workspace+'/'+"gg")
polygonBox=arcpy.MinimumBoundingGeometry_management(outputFeatureClass,'gg',"RECTANGLE_BY_AREA")
arcpy.AddField_management(outputFeatureClass,'angle',"DOUBLE")

#游标
cursor=arcpy.da.SearchCursor(polygonBox,['SHAPE@','ORIG_FID'])
cursorUpdate=arcpy.da.UpdateCursor(outputFeatureClass,['OID@','angle'])

#记录
recordDic={}
for row in cursor:
    points=row[0].getPart()[0]
    d1=distance(points[0],points[1])
    d2=distance(points[1],points[2])
    if d1>d2:
        recordDic[row[1]]=angle(points[0],points[1])
    else:
        recordDic[row[1]]=angle(points[1],points[2])
del cursor

#更新
for row in cursorUpdate:
    row[1]=recordDic[row[0]]
    cursorUpdate.updateRow(row)

转载自:https://blog.csdn.net/weixin_34397291/article/details/87412417

You may also like...