Viewshed Analysis

Be aware, for the following to work you need to have GRASS GIS installed on your system.

Load library, initialize the GRASS environment with default values (everything is stored as temporary files; you can change this if you want to come back to earlier stages of your work)

library(spgrass6)
loc <-  initGRASS("/usr/local/grass-7.0.0svn/", home=tempdir() ,mapset = "PERMANENT", override = TRUE)
execGRASS("g.proj", flags = c("c"), parameters = list(proj4="+init=epsg:32644"))
## Default region was updated to the new projection, but if you have multiple
## mapsets `g.region -d` should be run in each to update the region from the
## default
## Projection information updated
library(rgdal)
srtm <- readGDAL("./results/srtm.tif")
## ./results/srtm.tif has GDAL driver GTiff 
## and has 481 rows and 289 columns
writeRAST6(x = srtm, vname="dem", flags = c("overwrite"))
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## Raster map <dem> created.
## r.in.gdal complete.
## adjust Regions resolution!
execGRASS("g.region", parameters = list(raster = "dem",res = as.character(srtm@grid@cellsize[1])))
execGRASS("g.region", flags = c("p"))
## projection: 1 (UTM)
## zone:       44
## datum:      wgs84
## ellipsoid:  wgs84
## north:      1087615.18883952
## south:      654715.18883952
## west:       337168.57873155
## east:       597268.57873155
## nsres:      900
## ewres:      900
## rows:       481
## cols:       289
## cells:      139009

Now we are ready to perform our viewshed analysis.

# display the possible commands in grass
parseGRASS("r.viewshed")
## Command: r.viewshed 
## Description: Computes the viewshed of a point on an elevation raster map. 
## Keywords: Default format: NULL (invisible), vertical angle wrt viewpoint (visible). 
## Parameters:
##   name: input, type: string, required: yes, multiple: no
##   keydesc: name, keydesc_count: 1
## [Name of input elevation raster map]
##   name: output, type: string, required: yes, multiple: no
##   keydesc: name, keydesc_count: 1
## [Name for output raster map]
##   name: coordinates, type: float, required: yes, multiple: no
##   keydesc: east,north, keydesc_count: 2
## [Coordinates of viewing position]
##   name: observer_elevation, type: float, required: no, multiple: no
##   default: 1.75
##   keydesc: value, keydesc_count: 1
## [Viewing elevation above the ground]
##   name: target_elevation, type: float, required: no, multiple: no
##   default: 0.0
##   keydesc: value, keydesc_count: 1
## [Offset for target elevation above the ground]
##   name: max_distance, type: float, required: no, multiple: no
##   default: -1
##   keydesc: value, keydesc_count: 1
## [Maximum visibility radius. By default infinity (-1)]
##   name: refraction_coeff, type: float, required: no, multiple: no
##   default: 0.14286
## [Refraction coefficient]
##   name: memory, type: integer, required: no, multiple: no
##   default: 500
##   keydesc: value, keydesc_count: 1
## [Amount of memory to use in MB]
##   name: directory, type: string, required: no, multiple: no
## [Directory to hold temporary files (they can be large)]
## Flags:
##   name: c [Consider the curvature of the earth (current ellipsoid)] {FALSE}
##   name: r [Consider the effect of atmospheric refraction] {FALSE}
##   name: b [Output format is invisible = 0, visible = 1] {FALSE}
##   name: e [Output format is invisible = NULL, else current elev - viewpoint_elev] {FALSE}
##   name: overwrite [Allow output files to overwrite existing files] {FALSE}
##   name: help [Print usage summary] {FALSE}
##   name: verbose [Verbose module output] {FALSE}
##   name: quiet [Quiet module output] {FALSE}

Let us create viewshed maps for the peak data

peaks <- readOGR("./data/DSL250-Shp", "Peaks")
## OGR data source with driver: ESRI Shapefile 
## Source: "./data/DSL250-Shp", layer: "Peaks"
## with 21 features
## It has 6 fields
library(raster)
peaks <- spTransform(peaks, raster(srtm)@crs)

## for better orientation insert the boundary
boundary <- readOGR("./data/DSL250-Shp", "Boundary")
## OGR data source with driver: ESRI Shapefile 
## Source: "./data/DSL250-Shp", layer: "Boundary"
## with 504 features
## It has 5 fields
boundary <- spTransform(boundary, raster(srtm)@crs)
plot(boundary)
points(peaks, pch = 19, cex = .5)

## crop srtm to be relevant for dagoba data
# srtm_crop <- crop(raster(srtm), rbind(c(min(dagoba@bbox[1,]-10000),
#                                 max(dagoba@bbox[1,]+10000)),
#                               c(min(dagoba@bbox[2,]-10000),
#                                 max(dagoba@bbox[2,]+10000))
#                               ))      

srtm_crop <- raster(srtm)
writeRAST6(x = as(srtm_crop,"SpatialGridDataFrame"), vname="dem_crop", flags = c("overwrite"))
execGRASS("g.region", parameters = list(raster = "dem_crop"))

## viewshed for one point
co.peaks <- peaks@coords
execGRASS("r.viewshed", flags = c("overwrite","b"), parameters = list(input = "dem_crop",output = "view.peak",coordinates = co.peaks[1,]))
single.viewshed <- readRAST6("view.peak")
plot(raster(single.viewshed))
points(peaks[1,])

Question What might be wrong with this result?

execGRASS("r.viewshed", flags = c("overwrite","b"), parameters = list(input = "dem_crop",output = "view.peak",coordinates = co.peaks[1,], max_distance=50000))
single.viewshed <- readRAST6("view.peak")
plot(raster(single.viewshed))
points(peaks[1,],pch=19)

## viewshed for all points
## ------------------------------
## load basic raster from grass
dem_crop <- readRAST6(vname = "dem_crop")

## loop through all points, calculate viewshed, and write in a an raster brick object
cum.view.peaks <- brick(raster(dem_crop))

for (i in seq(1, length(co.peaks[,1]))) {
    execGRASS("r.viewshed"
              ,flags = c("overwrite","b")
              ,parameters = list(input = "dem_crop",
                                 output = "view.peak",
                                 coordinates = co.peaks[i,],
                                 max_distance=75000)
              )
    viewshed <- readRAST6("view.peak")
    cum.view.peaks[[i]] <- raster(viewshed)
    names(cum.view.peaks[[i]]) <- as.character(peaks$NAME[i])
    cat("iteration ", i, " of ", length(co.peaks[,1]),"\n")
}

plot(cum.view.peaks)

view.sum.peaks <- sum(cum.view.peaks)
plot(view.sum.peaks)

quantile(view.sum.peaks)
##   0%  25%  50%  75% 100% 
##   0    0    1   2  9 
##   1    2    3   4   5
rcl.peaks <- c(-Inf, quantile(view.sum.peaks)[2],1,
         quantile(view.sum.peaks)[2],quantile(view.sum.peaks)[3],2,
         quantile(view.sum.peaks)[3],quantile(view.sum.peaks)[4],3,
         quantile(view.sum.peaks)[4],quantile(view.sum.peaks)[5],4,
         quantile(view.sum.peaks)[5],+Inf,5
         )
rcl.peaks <- matrix(rcl.peaks, ncol = 3, byrow = TRUE)
view.sum.peaks2 <- reclassify(x = view.sum.peaks, rcl = rcl.peaks)
plot(view.sum.peaks2)
points(peaks, pch=19, cex=.5)

#library(rasterVis)
#levelplot(view.sum.peaks2)

Exercise

  1. There is something seriously wrong in this analysis. What is it?