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

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

MBR Tool

思路步骤
全部代码如下,作为个人记录,也分享给大家。。。
# -*- 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