CE-QUAL-W2 topographic shading

CE-QUALtopoShadingZenith.py

# Calculate 18 zenith angles (radians) for CE-QUAL-W2 Topographic Shading
# Use ArcGIS "Solar Radiation Graphics" Tool, output optional sky map
# Sky Size: 1440, Zenith Divisions: 360 (0.25 deg), Azimuth Divisions: 72 (5 deg)
# Tool returns a viewshed with as many bands as input points, and a sky map which
# has values starting at zero at the north outside edge, and spiraling clockwise inward
# However, the results are looking up at the sky so West/East are reversed
# Each sky map wedge is centered on the main direction
# Sky maps values divisible by 72 with no remainder (modulo or %), are from 357.5 to 2.5 azimuth
# Example: skymap value of 6210, where is it?
# Solution: 6210 % 72 = 18, There are 72 azimuth divisions so 18 is due west
#			6210 / 72 (integer division) = 86 so it's the 87th zenith loop
#			With 360 zenith divisions each division is 0.25 degrees so 21.75 degrees zenith
# This script executed in the arcpy window using results from the tool
# For simplicity, we will output all 72 zenith angles in their native reverse-order
# This will be post-processed in Excel, reversing the order and extracting every 4th result

import numpy as np
from scipy import ndimage as ndi

vsheds = arcpy.RasterToNumPyArray("vshed1440")
sky = arcpy.RasterToNumPyArray("skymap1440")
rad = np.deg2rad((sky/72 + 1)*0.25)
rad = rad * (rad > 0)
wedges = sky % 72
output = []
for v in vsheds:
	output.append(ndi.minimum(rad * v,wedges,range(72))*-1)
#np.savetxt('C:/Users/rsherrick/Desktop/zenith.csv',np.column_stack(output),delimiter=',')