I have a shapefile with >1200 polygons and wish to extract mean aspect (0-360) for all of them from an SRTM DEM (raster). Both are WGS84 GCS. I have found Python code to do this (here: https://geonet.esri.com/thread/47864), however get errors when trying to run.
Is there another way of doing this?
I have read other posts but have not fully understood what to do - I'm inexperienced in ArcGIS.
====================================================
SCRIPT:
CalcZonalMeanAspect.py
DESCRIPTION:
This script calculates the mean aspect azimuth (like 52 degrees) and the mean aspect
class (like 045 NE) for each zone (polygon) in a shapefile, or MDB or GDB feature class.
Mean aspect value fields are appended to the input's attribute table.
USER INPUTS:
> A polygon shapefile, or MDB or GDB feature class, containing the zones of interest
> A field name representing each polygon's unique ID
> An elevation raster
====================================================
Import modules and Spatial Analyst.
import arcpy, math, os
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
arcpy.AddMessage(" ")
arcpy.AddMessage(" ")
arcpy.AddMessage("Calc Zonal Mean Aspect was developed by Carl Beyerhelm, Coconino National Forest.")
arcpy.AddMessage(" ")
try:
# Accept user input.
inZone = arcpy.GetParameterAsText(0) ## The "zone" polygon shapefile or feature class
zoneField = arcpy.GetParameterAsText(1) ## The zone ID field
inElevation = arcpy.GetParameterAsText(2) ## The elevation raster
# Set the workspace to inZone's parent folder.
currentDir = inZone
scratchFolder = ""
while scratchFolder == "":
currentDir = os.path.dirname(currentDir)
if os.path.isdir(currentDir): ## If currentDir is a folder
if currentDir[-4:] == '.gdb': ## If currentDir is a file geodatabase
currentDir = os.path.dirname(currentDir) ## Retreat one more directory level
scratchFolder = currentDir
arcpy.env.workspace = scratchFolder ## Set the workspace to inZone's parent folder
# Set environment variables.
arcpy.env.mask = inZone
arcpy.env.extent = inZone
arcpy.env.snapRaster = inElevation
arcpy.env.overwriteOutput = True
# Generate a temporary aspect grid named asp from inElevation.
arcpy.AddMessage(" ")
arcpy.AddMessage("Generating an aspect grid...")
asp = arcpy.sa.Aspect(inElevation)
# asp.save("xxasp")
# del asp
# Set asp cells having a value of -1 (flat) to NULL.
arcpy.AddMessage(" ")
arcpy.AddMessage("Setting cell values with a FLAT aspect to NO DATA...")
aspNull = arcpy.sa.SetNull(asp, asp, "VALUE = -1")
# aspNull.save("xxaspnull")
# del aspNull
del asp
# Generate the cosine of aspect expressed in radians.
arcpy.AddMessage(" ")
arcpy.AddMessage("Getting the cosine of each cell's aspect expressed in radians...")
aspCos = arcpy.sa.Cos(aspNull * math.pi / 180.0)
# aspCos.save("xxaspcos")
# del aspCos
# Generate the sine of aspect expressed in radians.
arcpy.AddMessage(" ")
arcpy.AddMessage("Getting the sine of each cell's aspect expressed in radians...")
aspSin = arcpy.sa.Sin(aspNull * math.pi / 180.0)
# aspSin.save("xxaspsin")
# del aspSin
del aspNull
# Calculate the zonal sums of aspect cosine and aspect sine.
arcpy.AddMessage(" ")
arcpy.AddMessage("Calculating the zonal sums of aspect cosine and aspect sine...")
arcpy.sa.ZonalStatisticsAsTable(inZone, zoneField, aspCos, "xxSumCos.dbf", "DATA", "MEAN")
arcpy.sa.ZonalStatisticsAsTable(inZone, zoneField, aspSin, "xxSumSin.dbf", "DATA", "MEAN")
del aspCos
del aspSin
# Add fields to inZone's attribute table.
arcpy.AddMessage(" ")
arcpy.AddMessage("Adding fields Aspect_Azm and Aspect_Cls to " + os.path.basename(inZone) + "...")
arcpy.AddField_management(inZone, "MeanCos", "DOUBLE")
arcpy.AddField_management(inZone, "MeanSin", "DOUBLE")
arcpy.AddField_management(inZone, "Aspect_Azm", "DOUBLE")
arcpy.AddField_management(inZone, "Aspect_Cls", "TEXT", "", "", 6)
# Populate fields MeanCos and MeanSin with mean values from xxSumCos.dbf and xxSumSin.dbf.
arcpy.MakeTableView_management(inZone, "zoneView")
arcpy.MakeTableView_management("xxSumCos.dbf", "cosView")
arcpy.MakeTableView_management("xxSumSin.dbf", "sinView")
arcpy.AddJoin_management("zoneView", zoneField, "cosView", zoneField) ## Populate MeanCos values
arcpy.CalculateField_management("zoneView", "MeanCos", "[xxSumCos.Mean]")
arcpy.RemoveJoin_management("zoneView", "xxSumCos")
arcpy.AddJoin_management("zoneView", zoneField, "sinView", zoneField) ## Populate MeanSin values
arcpy.CalculateField_management("zoneView", "MeanSin", "[xxSumSin.Mean]")
arcpy.RemoveJoin_management("zoneView", "xxSumSin")
# Calculate the mean aspect azimuth for each zone.
arcpy.AddMessage(" ")
arcpy.AddMessage("Calculating the mean aspect azimuth for each zone in " + os.path.basename(inZone) + "...")
expression = "math.fmod(360 + (math.atan2(!MeanSin!, !MeanCos!)) * (180 / math.pi), 360)" ## This is the fundamental calculation
arcpy.CalculateField_management("zoneView", "Aspect_Azm", expression, "PYTHON")
# Calculate the mean aspect class for each zone.
arcpy.AddMessage(" ")
arcpy.AddMessage("Calculating the mean aspect class for each zone in " + os.path.basename(inZone) + "...")
arcpy.CalculateField_management("zoneView", "Aspect_Cls", '"' + '000' + ' NO' + '"') ## Calc everything to "000 NO" as a default
clsList = ["045 NE", "090 EA", "135 SE", "180 SO", "225 SW", "270 WE", "315 NW"] ## A list of aspect classes
low = 22.5
high = 67.5
fc = inZone
field = "Aspect_Azm"
delimField = arcpy.AddFieldDelimiters(fc, field) ## Construct a properly delimited field name
for cls in clsList:
expression = delimField + ' >= ' + str(low) + ' and ' + delimField + ' < ' + str(high) ## Select zones with an aspect azimuth between the
arcpy.SelectLayerByAttribute_management("zoneView", "NEW_SELECTION", expression) ## low and high bounds of the current aspect class
arcpy.CalculateField_management("zoneView", "Aspect_Cls", '"' + cls + '"', "PYTHON")
low = low + 45
high = high + 45
arcpy.AddMessage(" ")
arcpy.AddWarning("OK, done! Mean aspect azimuth and mean aspect class have been appended to:")
arcpy.AddWarning(" > " + inZone)
arcpy.AddMessage(" ")
except:
arcpy.AddError(" ")
arcpy.AddError("Something is cooked...")
arcpy.AddError(arcpy.GetMessages(2))
arcpy.AddError(" ")
finally:
# Clean up.
killFiles = ["xxSumCos.dbf", "xxSumSin.dbf"] ## Delete remnant .dbf files
for file in killFiles:
if arcpy.Exists(file):
arcpy.Delete_management(file)
killFields = ["MeanCos", "MeanSin"] ## Delete remnant fields
fieldList = arcpy.ListFields(inZone)
for field in fieldList:
if field.name in killFields:
arcpy.DeleteField_management(inZone, field.name)
killRasters = ["asp", "aspNull", "aspCos", "aspSin"] ## Dismiss remnant temporary rasters
for raster in killRasters:
if arcpy.Exists(raster):
del raster
arcpy.CheckInExtension("Spatial") ## Return SA license
arcpy.RefreshCatalog(scratchFolder) ## Refresh ArcCatalog view
