Great sample. Thought to add for reference if one wanted to calculate the output as Raster in Map Algebra

from arcpy.sa import *
import math
asp = Aspect(inElevation)
aspNull = SetNull(asp, asp, "VALUE = -1")
aspCos = Cos(aspNull * math.pi / 180.0)
aspSin = Sin(aspNull * math.pi / 180.0)
xxSumSin = ZonalStatistics(inZones,"OBJECTID", "aspSin", "MEAN", "DATA")
xxSumCos = ZonalStatistics(inZones,"OBJECTID", aspCos, "MEAN", "DATA")
asp_azm = (360+(ATan2(xxSumSin,xxSumCos)) * (180 / math.pi)) % 360.0