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:
Now you need to write
viewshed
andcan_see
. Theviewshed
function is going to be something involving the GRASS functions and it should, I guess, return a raster of 0 and 1:Then
can_view
takes two binary rasters with the same extent and resolution and returnsTRUE
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: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 simplemessage("Coord: ",i)
in the loop.