Line of Sight (LoS) Analysis: 3D Terrain Analysis (Part 2)

In my previous post on LoS Analysis, I have tried to explain briefly the basics of LoS in two dimensional space. Obviously real life problems are based on three dimensional terrains although basic concepts are all the same. In this second part I will try to adapt the same techniques with a few modifications for three dimensional terrains.

3D Terrain Visualization with R

One of the first differences in 3D LoS analysis is the terrain visualization. We can not use plot function for proper visualization is 3D. Fortunately R has all packages you need for any type of problem. I will use rgl package which can be downloaded using install.packages("rgl") command.

Once you have the rgl package, generating pseudo 3D terrains as we did for 2D is a trivial thing.

3D Terrain with RGL

You can use the following R script to generate your 3D terrains like above.

library(rgl)

# 3D Terrain Function
height <- function (x,y) {
  sin(x)+0.125*y*sin(2*x)+sin(y)+0.125*x*sin(2*y)+0.25
}

# Terrain boundaries -8<=x<=8 and -8<=y<=8
boundary <- c(-8,8)

# Terrain grid with a step size of 0.1 units
xy<-seq(from=boundary[1],to=boundary[2],by=0.1)

# Evaluate all heights for all grid points
z<-outer(xy,xy,height)

# A few visualization staff
zlim <- range(z)
zlen <- zlim[2] - zlim[1] + 1
colorlut <- terrain.colors(zlen) # height color lookup table
col <- colorlut[ z-zlim[1]+1 ] # assign colors to heights for each point

# Draw the terrain
rgl.open()
bg3d("gray")
rgl.surface(xy, xy, z,
            color=col)

A new function in this script is outer function which generates the product of a vector and a row-vector to have a matrix (product of a row-vector with a vector/column-vector is obviously a scalar value and named to be dot/inner product). The third parameter of the function provides us the mechanism to apply a given function (height in our case) for each element of this matrix. Obviously you can play with height function to have fancier 3D terrains and to have best visualization you may need viewpoint routine in rgl package .

LoS in 3D Terrain

Line of Sight analysis on 3D terrain uses the same principles as it does in 2D. Use the following R script to decide on status of a point (invisible, visible, visible but far away)

library(rgl)
##################
# Functions
##################
# 3D Terrain Function
height <- function (x,y) {
  sin(x)+0.125*y*sin(2*x)+sin(y)+0.125*x*sin(2*y)+0.25
}

# Linear Function
linear <- function (x, observer, target) {
  v <- observer - target

  y <- ((x - observer[1])/v[1])*v[2]+observer[2]
  z <- ((x - observer[1])/v[1])*v[3]+observer[3]

  data.frame(x=x,y=y, z=z)
}

# Linear Function
distance <- function (p0,p1) {
  sqrt(sum((p0-p1)^2))
}

##################
# Input
##################
# Observer location
observer<-c(10,10,1)

# Target on terrain
target <- c(5, 5, height(5,5))

# Max visible distance
maxVisibleDistance = 4

# Generate points with a step size of 0.1
x <- seq(from=min(observer[1],target[1]),
         to=max(observer[1],target[1]),
         by=0.1)

# All points on line
line <- linear(x, observer, target)

# Terrain Height
h <- height(line$x,line$y)

# LoS Analysis
aboveTerrain <- round((line$z-h),2) >= 0.1

# First Rule
visible <- !is.element(FALSE,aboveTerrain)
if (visible){
  # Second Rule
  d <- distance(observer, target)
  if(d <= maxVisibleDistance){
    status <- "LoS"
  }else{
    status <- "non-LoS due to Distance"
  }
}else{
  status <- "non-LoS due to Blocking"
}

Obviously there are a few changes in the script with compared to 2D version.  The first one is linear function(Code Lines 10-18). New version not only evaluates second (y) but also the third dimension (z). Notice that z is our height dimension by convention. We have also utilized data.frame function to concatenate all dimensions to form a table of point dimensions

The second difference is on height function (Code Lines 5-8). It is no longer a mapping from x to y but a mapping from x,y to z.

Rest of the 3D version of script is pretty much the same or trivial to discuss more.

Visualizing LoS on 3D Terrain

Until this point we have analyzed LoS of a single point on 2D-3D terrains. But usually network analists wish to know LoS map of the terrain with respect to a given observer.  In other words we need to visually understand which regions on 3D terrain are visible by the observer, invisible by the observer due to blocking, or further than the limit from the observer.

Here the LoS map of our pseudo 3D terrain with respect to an observer with a given set of coordinates and maximum service range(green vs yellow regions).

3D LoS Analysis

You can obtain this visualization using following R script.

library(rgl)
##################
# Functions
##################
# 3D Terrain Function
height <- function (point) {
  sin(point$x)+0.125*point$y*sin(2*point$x)+sin(point$y)+0.125*point$x*sin(2*point$y)+3
}

# Linear Function
linear <- function (px, observer, target) {
  v <- observer - target

  y <- ((px - observer[1])/v[1])*v[2]+observer[2]
  z <- ((px - observer[1])/v[1])*v[3]+observer[3]

  data.frame(x=px,y=y, z=z)
}

# Linear Function
distance <- function (terrain, observer) {
  sqrt((terrain$x-observer[1])^2+(terrain$y-observer[2])^2+(terrain$height-observer[3])^2)
}

LoS <- function(terrain, observer, maxVisibleDistance){
  status = c()
  for (i in seq(1:nrow(terrain))) {
    if (observer[1] == terrain$x[i] && observer[2] == terrain$y[i]){
      if(observer[3] >= terrain$height[i]){
        if (terrain$dist2observer[i] > maxVisibleDistance){
          status <- c(status,"yellow")
        }else{
          status <- c(status,"green")
        }
      }else{
        status <- c(status,"red")
      }
    }else{
      # All points on line
      line <- linear(seq(from=min(observer[1],terrain$x[i]),
                         to=max(observer[1],terrain$x[i]),
                         by=0.1),
                     observer,
                     c(terrain$x[i],terrain$y[i],terrain$height[i]))
      # Terrain Height
      h <- height(line)

      # LoS Analysis
      aboveTerrain <- round((line$z-h),2) >= 0.00

      visible <- !is.element(FALSE,aboveTerrain)
      if (visible){
        # Second Rule
        if(terrain$dist2observer[i] <= maxVisibleDistance){
          status <- c(status,"green")
        }else{
          status <- c(status,"yellow")
        }
      }else{
        status <- c(status,"red")
      }
    }
  }

  status
}

##################
# Input
##################
# Observer location
observer<-c(0.835597146302462, -1.71025141328573, 6)

# Max visible distance
maxVisibleDistance = 8

# Generate points with a step size of 0.1
x <- seq(from=-8,to=8,by=0.4)

xygrid <- expand.grid(x=x,
                      y=x)

terrain <- data.frame(xygrid,
                      height=height(xygrid)
                      )

terrain <- data.frame(terrain,
                      dist2observer=distance(terrain, observer)
                      )

terrain <- data.frame(terrain,
                      status = LoS(terrain, observer, maxVisibleDistance))

rgl.open()
rgl.surface(x, x,
            matrix(data=terrain$height,nrow=length(x),ncol=length(x)),
            col=matrix(data=terrain$status,nrow=length(x),ncol=length(x))
            )
bg3d("gray")
# Mark the observer
spheres3d(c(observer[1]),
          c(observer[3]),
          c(observer[2]),
          radius=0.5,
          color="white"
          )

rgl.viewpoint(-60,30)

For a better visualization R allows you to implement spinning 3D terrains using play3d function and record it in gif format using movie3d function as I did below.

3D LoS Analysis 360

Advertisement