Trigonometry to find landscape visibility in Python

Discussion created by dcn23 on Sep 11, 2013
Latest reply on Sep 12, 2013 by msayler
I am trying to find the degrees of visibility of natural environments from a number of viewpoints. I have a DEM, and a raster defining the natural environments. For each cell I have slope, height, aspect and distance from viewpoint. I want to calculate the angle between the viewpoint and each cell.

I originally did this by getting the x,y,z coordinates for the viewpoint and the both the bottom and top of a sloped cell. Then get the 3D distance between each coordinate and use these distances to get the angles. I would then sum for each cell to get total degrees of visibilty.

Problem is that this approach only works when the bearing between the viewpoint and the visible cell is in line with the x and y axis of the coordinate system. This because to calculate the points coordinates for the cell corners I am subtracting and adding the cell resolution/2 to the cell centroid.

Code below

def Cell_Loc (VPX, VPY, VPZ, cell_x, cell_y, cell_z, aspect, cell_resolution, bearing):
    import math
    from math import sqrt
    #Get change in height of cell using slope and cell resolution (trig)
    AspectTan = math.tan(math.radians(aspect))
    Opposite = (AspectTan * cell_res)

    #Get coordinates for cell corners
    CloseCornerX = cent_x - 2.5
    CloseCornerY = cent_y - 2.5
    FarCornerX = cent_x + 2.5
    FarCornerY = cent_y + 2.5
    CloseCornerZ = cent_z - (Opposite/2)
    FarCornerZ = cent_z + (Opposite/2)

#Get Distance between coordinates   
    VP = (VPX, VPY, VPZ) #data point 1
    cellcrner1 = (CloseCornerX, CloseCornerY, CloseCornerZ) #data point 2
    cellcrner2 = (FarCornerX, FarCornerY, FarCornerZ)    
    dist_to_far_corner = sqrt(sum( (VP - cellcrner2)**2 for VP, cellcrner2 in zip(VP, cellcrner2)))
    dist_to_close_corner = sqrt(sum( (VP - cellcrner1)**2 for VP, cellcrner1 in zip(VP, cellcrner1)))
    cell_dist = sqrt(sum( (cellcrner1 - cellcrner2)**2 for cellcrner1, cellcrner2 in zip(cellcrner1, cellcrner2)))

#Input above distances into ViewAngle to find angle 
def ViewAngle (a, b, c):
    "calculates the viewing angle of each visible cell"
    import math

    a2 = a*a
    b2 = b*b
    c2 = c*c
    VA = (a2-b2+c2)/ (2.0 * (a*c))
    rads = math.acos(VA)
    return math.degrees(rads)

How can I improve this code, or maybe another approach is a better way of doing it. The output I need is to say "for this viewpoint, there are X amount of degrees for visible natural environments". Again I have, DEM, visible areas from viewshed analysis, slope, aspect and distances between viewpoint and cell centroid. I use python programming language and have ArcGIS10.2 with arcpy