【ArcPy】空间相交求面积和
应用实例,求每宗地涉及开发边界的面积等。
实现思路(主干,非完整代码),谨供参考,不懂勿扰。
1 # coding=gbk 2 import arcpy 3 # import numpy 4 import sys 5 from arcpy import mapping as mapping 6 from arcpy import da as da 7 8 def main(): 9 target=arcpy.GetParameterAsText(0) 10 itsct=arcpy.GetParameterAsText(1) 11 12 tempItsct = arcpy.CreateScratchName("tempIntersectResult",data_type="FeatureClass") 13 inFcs=[target,itsct] 14 #获取图层Layer对象 15 targetIsLyr=isinstance(target, mapping.Layer) 16 targetLyr=target if targetIsLyr else mapping.Layer(target) 17 itsctIsLyr=isinstance(itsct, mapping.Layer) 18 itsctLyr=itsct if itsctIsLyr else mapping.Layer(itsct) 19 20 #获取图层线性单位,判断是不是米 21 targetLyrLinearUnit=arcpy.Describe(targetLyr).spatialReference.linearUnitName 22 itsctLyrLinearUnit=arcpy.Describe(itsctLyr).spatialReference.linearUnitName 23 if targetLyrLinearUnit!= 'Meter': 24 arcpy.AddError('>>>目标图层的线性单位不是米') 25 sys.exit() 26 if itsctLyrLinearUnit!= 'Meter': 27 arcpy.AddError('>>>相交图层的线性单位不是米') 28 sys.exit() 29 30 #获取图层名 31 targetLyrName=targetLyr.name 32 # itsctLyrName=targetLyr.name 33 #获取相交后的FID字段名,GDB字段名最大长度64 34 targetFldName=u"FID_{}".format(targetLyrName)[:65] 35 36 #执行相交操作 37 arcpy.Intersect_analysis(inFcs,tempItsct,'ONLY_FID','0.001','INPUT') 38 #转NumPy数组,使用Shape_Area 而非SHAPE@AREA,目的是保持与相交操作结果中字段Shape_Area的值一致。
39 arr= da.FeatureClassToNumPyArray(tempItsct,[targetFldName,'Shape_Area'])
40 #删除结果要素类 41 arcpy.Delete_management(tempItsct) 42 #分类求和 43 dic={} 44 for item in arr: 45 if item[0] in dic.keys(): 46 dic[item[0]]+=item[1] 47 else: 48 dic[item[0]]=item[1] 49 50 arcpy.AddMessage(dic) 51 52 53 if __name__ == '__main__': 54 main()