python - 三角函数寻找景观能见度

Question

我试图从多个角度找出自然环境的可见度。我有一个 DEM 和一个定义自然环境的栅格。对于每个单元格，我都有坡度、高度、方位角和距视点的距离。我想计算视点和每个单元格之间的角度。

我最初是通过获取视点的 x、y、z 坐标以及倾斜单元格的底部和顶部来完成此操作的。然后获取每个坐标之间的 3D 距离，并使用这些距离来获取角度。然后我会对每个单元格求和以获得总可见度。

问题是这种方法仅在视点和可见单元格之间的方位与坐标系的 x 轴和 y 轴一致时才有效。这是因为要计算单元格角点的坐标，我要减去单元格分辨率/2 并将其添加到单元格质心。

代码如下

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)   

我该如何改进这段代码，或者另一种方法可能是更好的方法。我需要的输出是说“对于这个观点，可见自然环境有 X 度数”。同样，我有 DEM、视域分析的可见区域、斜率、方面以及视点和单元质心之间的距离。我使用 python 编程语言并且有 ArcGIS10.2 和 arcpy

Answer 1

一个改进是使用坐标和内积而不是距离来计算角度:这将有助于避免一些(尽管不是全部)昂贵且困惑的平方根。所以如果你的观点是在v=(x,y,z) ，另外两个点是 P=(a,b,c)和 Q=(d,e,f) , 对角为 acos(sqrt(((V-P)(V-Q))^2/(V-P)^2(V-Q)^2))=acos(sqrt(<(V-P),(V-Q)>^2/<(V-P),(V-P)><(V-Q),(V-Q)>))其中(内部)产品是内部产品，即说 PQ = ad + be + cf .此外，一些产品和内部产品可以“回收”而不是重新计算。