我无法弄清楚如何在 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 几乎可以立即执行此操作,并且也可以更正确地执行此操作。
关于评估两个(纬度、经度)点之间的视线 (LOS) 的 R 代码,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21841387/