R – Calculating Interactions Between Two Viewsheds Using rgrass7

grass-gisrrgrass7viewshed

I have one shapefile with locations of an animal and one street layer.
First I created viewpoints on the street (every 10m) and created one cumulated tif ("street viewshed".
Now my plan is to create viewsheds from every single GPS point to calculate whether the "point viewshed" cuts the "street viewshed".
Not sure if ive explained it enough but the Output should be an R shapefile having the ID and the coordinates of the points and a third logical column (1=street viewshed gets cut by the point viewshed; 0=they don't interact).

I have started by using a for loop in R to create the viewsheds

view.brick <- brick(raster(dsm))
for (i in seq(1, length(coords[,1]))) {
  execGRASS("r.viewshed"
            ,flags = c("overwrite","b")
            ,parameters = list(input = "dsm",
                               output = "point_viewsheds",
                               coordinates = coords[i,]))
  viewshed <- readRAST("point_viewsheds")
  view.brick[[i]] <- raster(viewshed)
  cat("iteration ", i, " of ", length(coords[,1]),"\n")
}

With this I get it to create a .tmp file and folder for each point viewshed and make it create a rasterbrick. Unfortunately this brick is of the size of 3gb for only a 100 points (unfortunately I have over 100.000).

How do I achieve my target?

Best Answer

Here's an outline:

nc = nrow(coords)
overlap = rep(NA, nc) # make a vector to store the results

for(i in 1:nc){

  view_i = viewshed(dem, coords[i,])

  # optional: uncomment to save raster to an R RDS file:
  # saveRDS(view_i, paste0("view_",i,".rds"))

  overlap[i] = can_see(view_i, street_view)

}

Now you need to write viewshed and can_see. The viewshed function is going to be something involving the GRASS functions and it should, I guess, return a raster of 0 and 1:

viewshed = function(dem, pt){
  # grass stuff here
  #...
  #
  return(some_binary_raster)
}

Then can_view takes two binary rasters with the same extent and resolution and returns TRUE if there's any pixel with a 1 in any pair of cells. This is probably a one-liner that tests logical "AND" on the values from both rasters:

can_view = function(r1, r2){
   any(values(r1) & values(r2))
 }

Test this on a small subset of your coords before letting it loose on 100,000. You can use that to estimate how long it will take to do 100,000 coordinates. You could also add a progress bar (see ?txtProgressBar) or a simple message("Coord: ",i) in the loop.

Related Question