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]))) {
            ,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?

  • Does "every single GPS point" mean every point in the "animal" point data set?
    – Spacedman
    Commented Dec 14, 2021 at 16:36
  • Is the question essentially: Given a set of N points {a} and another set of M points {b}, compute the NxM rasters of the intersection of V(a_i) and V(b_j) where V(.) is the "visible from" function - ie the areas which can be seen by both points. How many points do you have in each set? How high resolution viewshed raster are you working with? How long does it take to do one calculation, for, say the first animal point and the first street point? Have you done that much? Could you show an example map?
    – Spacedman
    Commented Dec 14, 2021 at 16:41
  • Thanks for your reply! So on the one side i simply do have a viewshed generated from the street. It is a single raster file (10x10m) with values from 0-1 (0=areas are not visible, 1=very visible). On the other side ive got a shapefile that looks like ID geometry 1 x/y 2 x/y 3 x/y Now i want to create a third column indicating if the IDS are able to see the viewshed from the street (1=TRUE,0=FALSE). For that reason i have to create viewsheds for every id (with a smaller radius than for the street viewshed).
    – DennyCrane
    Commented Dec 14, 2021 at 17:28
  • @Spacedman With the code above it takes like 3 seconds to create a viewshed for one animal id. But the main problem is rather the size of the rasterbrick, my computer eventually wont handle it once i use thousands for points. With QGIS i am able to create the viewsheds around 30-40 times faster than with my current version of R Code. But QGIS also creates an endless amount of output files instead of a rasterbrick
    – DennyCrane
    Commented Dec 14, 2021 at 17:30
  • Why are you making such a brick? Why not make one viewshed single-layer raster for each point, intersect it with the street viewshed, and then its done with. Save it to a file if you want to keep it for later, you don't need it in R any more.
    – Spacedman
    Commented Dec 14, 2021 at 20:06

1 Answer 1


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

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.

  • thanks Spacedman. Worked perfectly fine the way you intended it to work! Unfortunately it takes more than 6 days to do it for a 100k points. I guess the main reason being that it creates a tmp file on my system for every viewshed. Tho it should still not take so long :D Also got another question: Ive got a dem, a ndsm (vegetation) and dsm (dem+ndsm). Obviously calculating viewsheds on a dsm isnt very fortunate since deer usually isnt standing on top of the trees. How do i create viewsheds of deer standing on the ground? My guess: dsm-ndsm for this coordinate, leaving the rest as it is?
    – DennyCrane
    Commented Dec 20, 2021 at 12:47
  • this will most likely not increase the calculation time of the 100k viewsheds tho :D Anyways, gotta try to find out how to make it faster
    – DennyCrane
    Commented Dec 20, 2021 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.