R code that evaluates line-of-sight (LOS) between two (lat, lon) points

Question

I'm having trouble figuring out how to calculate line-of-sight (LOS) between two (lat, lon) points, within R code. Any advice on how to approach this problem would be appreciated. I would like to use the R package - raster - for reading in the terrain elevation data. It seems the spgrass package could be leveraged (based on http://grass.osgeo.org/grass70/manuals/r.viewshed.html) but I wanted to avoid loading up a GIS. Thanks.

Answer 1

If you just want to know if point A can see point B then sample a large number of elevations from the line joining A to B to form a terrain profile and then see if the straight line from A to B intersects the polygon formed by that profile. If it doesn't, then A can see B. Coding that is fairly trivial. Conversely you could sample a number of points along the straight line from A to B and see if any of them have an elevation below the terrain elevation.

If you have a large number of points to compute, or if your raster is very detailed, or if you want to compute the entire area visible from a point, then that might take a while to run.

Also, unless your data is over a large part of the earth, convert to a regular metric grid (eg a UTM zone) and assume a flat earth.

I don't know of any existing package having this functionality, but using GRASS really isn't that much of a hassle.

Here's some code that uses raster and 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))])
}

Hopefully fairly self-explanatory but maybe needs some real documentation. I produced this with it:

which is a map of the points where a 50m high person can see a 2m high tower at the red dot. Yes, I got those numbers wrong when I ran it. It took about 20 mins to run on my 4 year old PC. I suspect GRASS could do this almost instantaneously and more correctly too.