setwd("YOUR-PATH")

raw_data <- read.table(file = "data/FS_UTM.csv",
                       header = TRUE,
                       sep = ";",
                       stringsAsFactors = FALSE,
                       quote = "")

## check for issues/problems
## str(raw_data)
## apply(X=raw_data[,4:20], MARGIN = 2, FUN=function(x){x[x>1]})
raw_data$settlement[raw_data$settlement > 1] <- 1

library(sp)
sites <- raw_data[(raw_data$settlement == 1 & raw_data$iron.age == 1),]
coordinates(sites) <- ~xUTM+yUTM
proj4string(sites) <- CRS("+init=epsg:32634")

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(rgrass7)

## Loading required package: XML

## GRASS GIS interface loaded with GRASS version: (GRASS not running)

loc <-  initGRASS(gisBase = "/usr/lib/grass72/", home=tempdir() ,mapset = "PERMANENT", override = TRUE)
execGRASS("g.proj", flags = c("c"), parameters = list(proj4="+init=epsg:32634"))

## 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 1705 rows and 3036 columns

library(raster)
srtm <- aggregate(raster(srtm), fact = 10)
library(maptools)
srtm <- as(object = srtm, Class = "SpatialGridDataFrame")

writeRAST(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%

Now, we need to adjust the region's 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:       34
## datum:      wgs84
## ellipsoid:  wgs84
## north:      5155478.7961
## south:      5001578.7961
## west:       433121.3283
## east:       706721.3283
## nsres:      900
## ewres:      900
## rows:       171
## cols:       304
## cells:      51984

Now we are ready to perform our viewshed analysis. First, display the possible commands of the respetive GRASS GIS module

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}

execGRASS("r.viewshed", flags = c("overwrite","b"),
          parameters = list(input = "dem",
                            output = "viewshed.single",
                            coordinates = sites@coords[1,]
          )
)

## Computing events...
##    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%
## Computing visibility...
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
## Writing output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%

single.viewshed <- readRAST("viewshed.single")

## Creating BIL support files...
## Exporting raster as integer values (bytes=4)
##    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%

execGRASS("r.viewshed", flags = c("overwrite","b"),
          parameters = list(input = "dem",
                            output = "viewshed.single",
                            coordinates = sites@coords[1,],
                            max_distance=50000
                            
          )
)

## Computing events...
##    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%
## Computing visibility...
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
## Writing output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%

single.viewshed <- readRAST("viewshed.single")

## Creating BIL support files...
## Exporting raster as integer values (bytes=4)
##    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%

Again a loop.

Procedure:

• for all points
• take i-th coordinate
• calculate viewshed
• save result in raster brick
• optional: plot progress information

view <- srtm
view@data$band1 <- 0

writeRAST(x = view,
          vname="view",
          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%

for (i in seq(1, length(sites@coords[,1]))) {
    execGRASS("r.viewshed",
              flags = c("overwrite","b", "quiet"),
              parameters = list(input = "dem",
                                output = "view_tmp",
                                coordinates = sites@coords[i,],
                                max_distance = 40000,
                                memory = 1000)
              )
    execGRASS("r.mapcalc",
              parameters = list(expression = paste0("'view' = 'view' + view_tmp")
                                ),
              flags = c("overwrite", "quiet")
              )
    #cat("iteration ", i, " of ", length(sites@coords[,1]),"\n")
}

view <- readRAST("view")

## Creating BIL support files...
## Exporting raster as floating values (bytes=8)
##    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%

plot(raster(view))
points(sites)

The Euclidean distance between points i and j is the length of the line segment connecting them.

\[ \begin{aligned} d(i,j) = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2} \end{aligned} \]

Manual calculation:

sqrt((sites@coords[1,1]-sites@coords[2,1])^2 + (sites@coords[1,2]-sites@coords[2,2])^2)

## [1] 41668.6

The Euclidean distance between points i and j is the length of the line segment connecting them.

\[ \begin{aligned} d(i,j) = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2} \end{aligned} \]

With own function:

e.d <- function(x, y) {
    sqrt(sum((x[1] - y[1])^2 + (x[2] - y[2])^2))
}
e.d(sites@coords[1,],sites@coords[2,])

## [1] 41668.6

The Euclidean distance between points i and j is the length of the line segment connecting them.

\[ \begin{aligned} d(i,j) = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2} \end{aligned} \]

With own function - for all connections:

x_comb <- combn(x = sites@coords[,1],
                m = 2)
y_comb <- combn(x = sites@coords[,2],
                m = 2)

ed_sites <- c()
for (i in 1:length(x_comb)/2) {
    ed_sites[i] <- e.d(cbind(x_comb[1,][i],y_comb[1,][i]), cbind(x_comb[2,][i],y_comb[2,][i]))
    }
ed_sites[1]

## [1] 41668.6

The Euclidean distance between points i and j is the length of the line segment connecting them.

\[ \begin{aligned} d(i,j) = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2} \end{aligned} \]

Using packages:

library(fields)
ed_sites2 <- rdist(x1 = sites@coords,
                   x2 = sites@coords
                   )
ed_sites2[2]

## [1] 41668.6

library(rgeos)
ed_sites3 <- gDistance(sites, byid = TRUE)
ed_sites3[2]

## [1] 41668.6

First: define cost functions

tobler1993a <- function(s){6 * exp(-3.5 * abs(s + 0.05))}      # km/h
tobler1993b <- function(s){0.36 * exp(-3.5 * abs(s + 0.05))}   # m/min

We use the gdistance package (link + nice vignette)

1. Auxilliary function (difference vector = slope)
2. Transitional object (which cells are connected in which direction (4,8,16)?)
3. Geocorrection (because a diagonal path is longer than a horizontal/vertical path)
4. Adjacency matrix (check object, i.e. use only those cells that are adjacent; since matrix values for non-adjacent cells are 0 and division by zero leads to Inf)
5. Apply cost function on adjacent cells (cell value = the difference in elevation from one to the other)

check this: transform cost to conductivity; conductivity=cost/dist; time=1/conductivity; we need the geocorection twice because

library(raster)
srtm <- raster(x = "results/srtm.tif")
srtm <- aggregate(srtm, fact = 10)

library(gdistance) # 
hdiff <- function(x){(x[2]-x[1])}

hd <- transition(srtm,hdiff,8,symm=FALSE)

## Warning in .TfromR(x, transitionFunction, directions, symm): transition
## function gives negative values

slope <- geoCorrection(hd,scl=FALSE)
adj <- adjacent(x=srtm, cells=1:ncell(srtm), direction=8)
cost <- slope       
cost[adj] <- tobler1993a(slope[adj])
conduct <- geoCorrection(cost, scl=FALSE)                                        

It's about time to calculate the shortest path

tmp1 <- shortestPath(conduct,
                     origin = sites@coords[1,],
                     goal = sites@coords[2,],
                     output="SpatialLines"
                     )

tmp2 <- shortestPath(conduct,
                     origin = sites@coords[2,],
                     goal = sites@coords[1,],
                     output="SpatialLines"
                     )

and plot

plot(srtm,
     main = "Least-cost path\nbased on Tobler's function for walking speed"
     )

lines(tmp1, col = "red")
lines(tmp2, col = "blue")
points(sites@coords[1,], pch = 19, cex = .5)
points(sites@coords[2,], pch = 19, cex = .5)
text(sites@coords[1,1], sites@coords[1,2] - 2500, "1st")
text(sites@coords[2,1], sites@coords[2,2] + 2500, "2nd")

legend("bottom",
       c("1st - 2nd", "2nd - 1st"),
       lty = c(1,1),
       col = c("red", "blue", "black"),
       bty = "n"
       )

Let us have a look at the elevation profiles of the calculated paths and compare them to the euclidean path

First, we need to create a spatial line between the 1st and 2nd point

tmp1_euc <- SpatialLines(
    list(
        Lines(
            Line(
                coords = (rbind(sites@coords[1,],
                                sites@coords[2,]
                                )
                )
            ),
            ID="1")
    ),
    proj4string = srtm@crs
)

Let us have a look at the elevation profiles of the calculated paths and compare them to the euclidean path

Now, we can extract the raster values along the spatial lines

p_tmp1_euc <- extract(x=srtm,
                      y=tmp1_euc, along=TRUE
                      )

p_tmp1 <- extract(x=srtm,
                  y=tmp1, along=TRUE
                  )
p_tmp2 <- extract(x=srtm,
                  y=tmp2, along=TRUE
                  )

Let us have a look at the elevation profiles of the calculated paths and compare them to the euclidean path

And now the plot

par(mfrow = c(3,1))
plot(p_tmp1_euc[[1]], type="l",
     main = "Elevation profile \nfrom 1st to 2nd site (euclidean)"
     )
plot(p_tmp1[[1]], type="l",
     main = "Elevation profile \nfrom 1st to 2nd site (least cost)"
     )
plot(p_tmp2[[1]], type="l",
     main = "Elevation profile \nfrom 2nd to 1st site (least cost)"
     )

Let us have a look at the elevation profiles of the calculated paths and compare them to the euclidean path

And now the plot

First we need a combination object that calculates all the possible connections (and optionally a name object in order to identify the paths afterwards)

x <- sites@coords[,1][1:20] # at the moment only for 20 objects
y <- sites@coords[,2][1:20] # at the moment only for 20 objects

x_comb <- combn(x,2)
y_comb <- combn(y,2)

x_name <- data.frame(sites@data$Nr.[match(x_comb,sites@coords[,1])])

This will be a rather long loop…

• for all points
• get odd and even numbers for directions, i.e. from odd to even and from even to odd
• select the odd/even coordinates from the combination object and store in variable
• calculate path from odd/even to even/odd
• write an ID to be able to identify the resulting lines
• in the first run of the loop: create the spatial data frame (sdf), i.e. every row is a calculated line
• in the subsequent runs of the loop: append the new lines to the sdf
• in the first run of the loop: create a dataframe (df) with columns "ID", "FROM", "TO"
• in the subsequent runs of the loop: append the new data

for(i in 1:(length(x_comb)/2)){
    i1 <- i*2-1 # odd numbers
    i2 <- i*2 # even numbers
    s <- c(x_comb[i1],y_comb[i1])
    z <- c(x_comb[i2],y_comb[i2])
    sz <- shortestPath(conduct, s, z, output="SpatialLines") # calculate the shortest path
    zs <- shortestPath(conduct, z, s, output="SpatialLines") # calculate the shortest path
    sz@lines[[1]]@ID <- as.character(paste(x_name[i1,1],x_name[i2,1]))
    zs@lines[[1]]@ID <- as.character(paste(x_name[i2,1],x_name[i1,1]))
    
    if(i==1){sdf <-rbind(sz,zs)}
    if(i>1){sdf <- rbind(sdf,sz,zs,
                         makeUniqueIDs = TRUE)
    }

    if(i==1){df <- cbind(ID = c(1,2),
                         FROM = c(paste(x_name[i1,1]),paste(x_name[i2,1])),
                         TO = c(paste(x_name[i2,1]),paste(x_name[i1,1]))
                         )
    }
    if(i>1){df <- cbind(ID = c(df[,1],i1,i2),
                        FROM = c(df[,2],paste(x_name[i1,1]),paste(x_name[i2,1])),
                        TO = c(df[,3],paste(x_name[i2,1]),paste(x_name[i1,1]))
                        )
    }
}

Create a nice, all integrating SpatialLinesDataFrame object, calculate the length of the lines, and write a shapefile

lcp_df <- as.data.frame(df)
# colnames(lcp_df) <- c("ID","START","TARGET")

lcp_sldf <- SpatialLinesDataFrame(sdf,
                                  lcp_df,
                                  match.ID = FALSE
                                  )

lcp_sldf@data$DIST <- SpatialLinesLengths(lcp_sldf)

writeOGR(lcp_sldf,
         "results/",
         "Least_cost_paths",
         driver = "ESRI Shapefile",
         overwrite_layer = TRUE
         )