arcgis python添加几何属性

import arcpy
import numpy
import mathdef AddGeometryAttributes(fc, geomProperties, lUnit, aUnit, cs):
"""-------------------------------------------------------------------------
Tool: Add Geometry Attributes (Data Management Tools)
Source Name: AddGeometryAttributes.py
Version: ArcGIS 10.2.1
Author: Esri, Inc.
Usage: arcpy.AddGeometryAttributes_management(
Input_Features,
Geometry_Properties,
{Length_Unit},
{Area_Unit},
{Coordinate_System})
Required Arguments: Input Features
Geometry Properties
Optional Arguments: Length Unit
Area Unit
Coordinate System
Description: Adds attribute fields to the input features containing
measurements and coordinate properties of the feature
geometries (for example, length or area).
Updated: Not yet.
------------------------------------------------------------------------""" desc = arcpy.Describe(fc)
hasZ, hasM = desc.hasZ, desc.hasM
if cs:
sr = arcpy.SpatialReference()
sr.loadFromString(cs)
try:
srMetersPerUnit = sr.metersPerUnit
except:
srMetersPerUnit = 1
else:
try:
srMetersPerUnit = desc.spatialReference.metersPerUnit
except:
srMetersPerUnit = 1
shapeDict = {"POINT":1,
"MULTIPOINT":1.5,
"POLYLINE":2,
"POLYGON":3,}
shapeDim = shapeDict[str(desc.shapeType).upper()]
del desc fields = CreateOutputFields(fc, geomProperties, hasZ, hasM) arcpy.SetProgressor("STEP", arcpy.GetIDMessage(86174), 0,
int(arcpy.GetCount_management(fc).getOutput(0)), 1) hasNulls = False
# Calculate geometry properties into new fields
with arcpy.da.UpdateCursor(fc,fields + ["SHAPE@"],"",cs) as ucur:
ucurFields = ucur.fields
global ucurFields
for row in ucur:
global row
geom = row[ucurFields.index("SHAPE@")]
if geom:
if shapeDim == 1:
if "POINT_X_Y_Z_M" in geomProperties:
Update("POINT_X", geom.firstPoint.X)
Update("POINT_Y", geom.firstPoint.Y)
if hasZ:
Update("POINT_Z", geom.firstPoint.Z)
if hasM:
Update("POINT_M", geom.firstPoint.M)
if shapeDim>1:
if "PART_COUNT" in geomProperties:
Update("PART_COUNT", geom.partCount)
if "CENTROID" in geomProperties:
Update("CENTROID_X", geom.trueCentroid.X)
Update("CENTROID_Y", geom.trueCentroid.Y)
if hasZ:
Update("CENTROID_Z", geom.trueCentroid.Z)
if hasM:
Update("CENTROID_M", geom.trueCentroid.M)
if "EXTENT" in geomProperties:
Update("EXT_MIN_X", geom.extent.XMin)
Update("EXT_MIN_Y", geom.extent.YMin)
Update("EXT_MAX_X", geom.extent.XMax)
Update("EXT_MAX_Y", geom.extent.YMax)
if shapeDim>=2:
if "POINT_COUNT" in geomProperties:
Update("PNT_COUNT", geom.pointCount)
if "LINE_START_MID_END" in geomProperties:
Update("START_X", geom.firstPoint.X)
Update("START_Y", geom.firstPoint.Y)
if shapeDim ==2:
midPoint = geom.positionAlongLine(0.5,
True).firstPoint
else:
line = arcpy.Polyline(geom.getPart(0), "#",
hasZ, hasM)
if line.length > 0:
midPoint = line.positionAlongLine(0.5,
True).firstPoint
else:
hasNulls = True
del line
Update("MID_X", midPoint.X)
Update("MID_Y", midPoint.Y)
Update("END_X", geom.lastPoint.X)
Update("END_Y", geom.lastPoint.Y)
if hasZ:
Update("START_Z", geom.firstPoint.Z)
Update("MID_Z", midPoint.Z)
Update("END_Z", geom.lastPoint.Z)
if hasM:
Update("START_M", geom.firstPoint.M)
Update("MID_M", midPoint.M)
Update("END_M", geom.lastPoint.M)
del midPoint
if "CENTROID_INSIDE" in geomProperties:
Update("INSIDE_X", geom.centroid.X)
Update("INSIDE_Y", geom.centroid.Y)
if hasZ:
Update("INSIDE_Z", geom.centroid.Z)
if hasM:
Update("INSIDE_M", geom.centroid.M)
if shapeDim==2:
if "LINE_BEARING" in geomProperties:
lat1 = geom.firstPoint.Y
lon1 = geom.firstPoint.X
lat2 = geom.lastPoint.Y
lon2 = geom.lastPoint.X
Update("BEARING", (90 - math.degrees(math.atan2(lat2-lat1, lon2-lon1))) % 360)
del lat1, lon1, lat2, lon2
if "LENGTH" in geomProperties:
Update("LENGTH", ConvertFromMeters("LINEAR",
geom.length,
lUnit,
srMetersPerUnit))
if "LENGTH_3D" in geomProperties:
Update("LENGTH_3D", ConvertFromMeters("LINEAR",
geom.length3D,
lUnit,
srMetersPerUnit))
if "LENGTH_GEODESIC" in geomProperties:
Update("LENGTH_GEO", ConvertFromMeters("LINEAR",
geom.getLength("PRESERVE_SHAPE"),
lUnit,
srMetersPerUnit))
if shapeDim==3:
if "PERIMETER_LENGTH" in geomProperties:
Update("PERIMETER", ConvertFromMeters("LINEAR",
geom.length,
lUnit,
srMetersPerUnit))
if "AREA" in geomProperties:
Update("POLY_AREA", ConvertFromMeters("AREA",
geom.area,
aUnit,
srMetersPerUnit))
if "AREA_GEODESIC" in geomProperties:
Update("AREA_GEO", ConvertFromMeters("AREA",
geom.getArea("PRESERVE_SHAPE"),
aUnit,
srMetersPerUnit))
if "PERIMETER_LENGTH_GEODESIC" in geomProperties:
Update("PERIM_GEO", ConvertFromMeters("LINEAR",
geom.getLength("PRESERVE_SHAPE"),
lUnit,
srMetersPerUnit))
ucur.updateRow(row)
else:
hasNulls = True
arcpy.SetProgressorPosition()
if hasNulls:
arcpy.AddIDMessage("WARNING", 957)def CreateOutputFields(fc, geomProperties, hasZ, hasM):
propDict = {"POINT_X_Y_Z_M": ["POINT_X",
"POINT_Y",
"POINT_Z",
"POINT_M"],
"PART_COUNT": ["PART_COUNT"],
"CENTROID": ["CENTROID_X",
"CENTROID_Y",
"CENTROID_Z",
"CENTROID_M"],
"EXTENT": ["EXT_MIN_X",
"EXT_MIN_Y",
"EXT_MAX_X",
"EXT_MAX_Y"],
"POINT_COUNT": ["PNT_COUNT"],
"LINE_START_MID_END": ["START_X",
"START_Y",
"START_Z",
"START_M",
"MID_X",
"MID_Y",
"MID_Z",
"MID_M",
"END_X",
"END_Y",
"END_Z",
"END_M"],
"LINE_BEARING": ["BEARING"],
"CENTROID_INSIDE": ["INSIDE_X",
"INSIDE_Y",
"INSIDE_Z",
"INSIDE_M"],
"LENGTH": ["LENGTH"],
"PERIMETER_LENGTH": ["PERIMETER"],
"AREA": ["POLY_AREA"],
"LENGTH_GEODESIC": ["LENGTH_GEO"],
"AREA_GEODESIC": ["AREA_GEO"],
"LENGTH_3D": ["LENGTH_3D"],
"PERIMETER_LENGTH_GEODESIC":["PERIM_GEO"],
}
if not hasZ:
propDict["POINT_X_Y_Z_M"].remove("POINT_Z")
propDict["CENTROID"].remove("CENTROID_Z")
propDict["CENTROID_INSIDE"].remove("INSIDE_Z")
propDict["LINE_START_MID_END"].remove("START_Z")
propDict["LINE_START_MID_END"].remove("MID_Z")
propDict["LINE_START_MID_END"].remove("END_Z")
if not hasM:
propDict["POINT_X_Y_Z_M"].remove("POINT_M")
propDict["CENTROID"].remove("CENTROID_M")
propDict["CENTROID_INSIDE"].remove("INSIDE_M")
propDict["LINE_START_MID_END"].remove("START_M")
propDict["LINE_START_MID_END"].remove("MID_M")
propDict["LINE_START_MID_END"].remove("END_M") addList = []
geomPropertiesList = []
currentFields = [field.name for field in arcpy.ListFields(fc)]
for prop in geomProperties:
for field in propDict[prop.upper()]:
geomPropertiesList.append(field)
if field not in currentFields:
addList.append(field)
else:
arcpy.AddIDMessage("WARNING", 1097, field) if addList:
narray = numpy.array([], numpy.dtype([("_ID", numpy.int)] +
[(n, numpy.float) for n in addList]))
arcpy.da.ExtendTable(fc, "OID@", narray, "_ID") return geomPropertiesListdef ConvertFromMeters(type, value, unit, metersPerUnit):
if not unit:
return value
else:
distanceUnitInfo = {"METERS": 1.0,
"FEET_US": 0.304800609601219,
"NAUTICAL_MILES": 1852.0,
"MILES_US": 1609.34721869444,
"KILOMETERS": 1000.0,
"YARDS": 0.914401828803658,} areaUnitInfo = {"ACRES": 4046.86,
"HECTARES": 10000.0,
"SQUARE_METERS": 1.0,
"SQUARE_FEET_US": 0.09290341161327473,
"SQUARE_NAUTICAL_MILES": 3429904.0,
"SQUARE_MILES_US": 2589998.4703195295,
"SQUARE_KILOMETERS": 1000000.0,
"SQUARE_YARDS": 0.8361307045194741,}
if type == "LINEAR":
return (value*metersPerUnit)/distanceUnitInfo[unit]
if type == "AREA":
return (value*math.pow(metersPerUnit,2))/areaUnitInfo[unit]def Update(field, value):
if value:
if not math.isnan(value):
row[ucurFields.index(field)] = value#run the script
if __name__ == '__main__':
# Get Parameters
fc = arcpy.GetParameterAsText(0)
if arcpy.GetParameterAsText(1).find(";") > -1:
geomProperties = arcpy.GetParameterAsText(1).upper().split(";")
else:
geomProperties = [arcpy.GetParameterAsText(1).upper()]
lUnit = arcpy.GetParameterAsText(2)
aUnit = arcpy.GetParameterAsText(3)
cs = arcpy.GetParameterAsText(4)
if not cs:
cs = arcpy.env.outputCoordinateSystem # Run the main script
AddGeometryAttributes(fc, geomProperties, lUnit, aUnit, cs)