ArcGIS mean aspect of polygons from SRTM

Question

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

What precisely have you tried? If you include the code you ran and the errors it gave then we may have a starting point to help. — PolyGeo, Dec 28 '16 at 19:43
@PolyGeo - I have tried to run the python script attached to the link above and get the following error: Something is cooked... ERROR 010328: Syntax error at or near symbol (. ERROR 010267: Syntax error in parsing grid expression. Failed to execute (Aspect). I'm using a shapefile of distinct polygons, and an SRTM elevation model — Darren J, Dec 29 '16 at 11:47
Requiring potential answerers to follow a link to find code and read comments to extract error messages is nowhere near as helpful as including a code snippet that illustrates where you are stuck and any errors produced in your question. — PolyGeo, Dec 29 '16 at 12:00
@PolyGeo - My apologies. I have attached the original code. I have been running it within the ArcMap shell, which does not (to my knowledge) highlight where errors occur. Added to 'GIS' topic as involves shapefiles etc, and is run within ArcMap using a built-in toolbox. Hope this is better? — Darren J, Dec 29 '16 at 12:54
I'll re-open for now but usually help is obtained only by extracting a code snippet rather than presenting code in its entirety so you may still struggle to attract potential answerers. — PolyGeo, Dec 29 '16 at 13:06

score 1 · Answer 1 · edited Apr 13 '17 at 12:33

1

A suitable answer may be present in Is there a zonal statistics tool to calculate median and mode raster values for each polygon?

It mentions median and mode in heading but the answers include mean calculation.

edited Apr 13 '17 at 12:33

Community

1

answered Dec 28 '16 at 17:09

lpreston

55
5

ArcGIS mean aspect of polygons from SRTM

====================================================

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.

Accept user input.

Set the workspace to inZone's parent folder.

Set environment variables.

Generate a temporary aspect grid named asp from inElevation.

asp.save("xxasp")

del asp

Set asp cells having a value of -1 (flat) to NULL.

aspNull.save("xxaspnull")

del aspNull

Generate the cosine of aspect expressed in radians.

aspCos.save("xxaspcos")

del aspCos

Generate the sine of aspect expressed in radians.

aspSin.save("xxaspsin")

del aspSin

Calculate the zonal sums of aspect cosine and aspect sine.

Add fields to inZone's attribute table.

Populate fields MeanCos and MeanSin with mean values from xxSumCos.dbf and xxSumSin.dbf.

Calculate the mean aspect azimuth for each zone.

Calculate the mean aspect class for each zone.

Clean up.

1 Answers1