r - R代码，用于评估两个(纬度，经度)点之间的视线(LOS)

我在弄清楚如何计算R代码中两个(纬度，经度)点之间的视线(LOS)时遇到了麻烦。任何有关如何解决此问题的建议，将不胜感激。我想使用R包-栅格-读取地形高程数据。似乎可以利用spgrass软件包(基于http://grass.osgeo.org/grass70/manuals/r.viewshed.html)，但我想避免加载GIS。谢谢。

最佳答案

如果只想知道点A是否可以看到点B，则从连接A到B的直线上采样大量高程以形成地形轮廓，然后查看从A到B的直线是否与该轮廓所形成的多边形相交。如果不是，则A可以看到B。这是相当琐碎的编码。相反，您可以沿A到B的直线对许多点进行采样，看是否有一个点的海拔高度低于地形高度。

如果要计算的点很多，或者栅格非常详细，或者要计算从一个点可见的整个区域，则可能需要一段时间才能运行。

另外，除非您的数据覆盖地球的大部分，否则请转换为常规的公制网格(例如，UTM区域)，并假设地球是平坦的。

我不知道任何现有的具有此功能的软件包，但是使用GRASS确实没有那么麻烦。

这是一些使用raster和plyr的代码:

cansee <- function(r, xy1, xy2, h1=0, h2=0){
### can xy1 see xy2 on DEM r?
### r is a DEM in same x,y, z units
### xy1 and xy2 are 2-length vectors of x,y coords
### h1 and h2 are extra height offsets
###  (eg top of mast, observer on a ladder etc)
    xyz = rasterprofile(r, xy1, xy2)
    np = nrow(xyz)-1
    h1 = xyz$z[1] + h1
    h2 = xyz$z[np] + h2
    hpath = h1 + (0:np)*(h2-h1)/np
    return(!any(hpath < xyz$z))
}

viewTo <- function(r, xy, xy2, h1=0, h2=0, progress="none"){
    ## xy2 is a matrix of x,y coords (not a data frame)
    require(plyr)
    aaply(xy2, 1, function(d){cansee(r,xy,d,h1,h2)}, .progress=progress)
}

rasterprofile <- function(r, xy1, xy2){
### sample a raster along a straight line between two points
### try to match the sampling size to the raster resolution
    dx = sqrt( (xy1[1]-xy2[1])^2 + (xy1[2]-xy2[2])^2 )
    nsteps = 1 + round(dx/ min(res(r)))
    xc = xy1[1] + (0:nsteps) * (xy2[1]-xy1[1])/nsteps
    yc = xy1[2] + (0:nsteps) * (xy2[2]-xy1[2])/nsteps
    data.frame(x=xc, y=yc, z=r[cellFromXY(r,cbind(xc,yc))])
}

希望很不言自明，但可能需要一些真实的文档。我用它制作了这个:

这是一个50m高的人可以在红点看到2m高的塔的点的 map 。是的，我在运行时弄错了这些数字。在我4岁的PC上运行大约需要20分钟。我怀疑GRASS几乎可以立即且更正确地执行此操作。