What is the shortest path from Colombo to Kandy?

library(rgdal)
cities <- readOGR(dsn = "./data/DSL250-Shp", layer = "Cities")
## OGR data source with driver: ESRI Shapefile 
## Source: "./data/DSL250-Shp", layer: "Cities"
## with 27 features
## It has 6 fields
str(cities@data)
## 'data.frame':    27 obs. of  6 variables:
##  $ AREA     : num  13432220 8275633 2041111 4956923 23486170 ...
##  $ PERIMETER: num  17189 12217 5707 12061 21727 ...
##  $ CITIES_  : num  2 3 4 5 6 7 8 9 10 11 ...
##  $ CITIES_ID: num  2 12 14 17 20 24 29 30 31 33 ...
##  $ CODE     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ NAME     : Factor w/ 25 levels "AMPARA","ANURADHAPURA",..: 9 13 19 15 25 2 21 4 4 4 ...
## as polygons
colombo <- cities[cities@data$NAME=="COLOMBO",]
kandy <- cities[cities@data$NAME=="KANDY",]

## as points
library(rgeos)
p.colombo <- gCentroid(cities[cities@data$NAME=="COLOMBO",])
p.kandy <- gCentroid(cities[cities@data$NAME=="KANDY",])
 
par(mfrow = c(1,2))
plot(colombo, main = "Colombo")
points(p.colombo)
plot(kandy, main = "Kandy")
points(p.kandy)

## fancy
## install.packages("plotGoogleMaps")
## library(plotGoogleMaps)
## plotGoogleMaps(p.kandy)

Euclidean Distance

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} \]

## Euclidean Distance - manual
sqrt((p.kandy@coords[1]-p.colombo@coords[1])^2+(p.kandy@coords[2]-p.colombo@coords[2])^2)
## [1] 93816.32
## Euclidean Distance - with own function
e.d <- function(x, y) {
    sqrt(sum((x - y)^2))
}
e.d(p.kandy@coords,p.colombo@coords)
## [1] 93816.32
## using package
library(fields)
rdist(p.kandy@coords,p.colombo@coords)
##          [,1]
## [1,] 93816.32
library(rgeos)
gDistance(p.kandy,p.colombo)
## [1] 93816.32

Least Cost Distance

## load srtm data
library(raster)
srtm <- raster(x = "./results/srtm.tif")

### LEAST COST PATH ###
#######################

nz <- 8   # neigbourhood number: 4,8,16
ras <- focal(srtm, w=matrix(1/9,nrow=3,ncol=3), NAonly=TRUE) # define the moving window and store it in the raster object; it is not allowed to chose the moving window too large because than a "jump" over high cost/ low conductivity cells would be possible
plot(ras, col = gray.colors(25, start = 0.97, end = 0.4))

# cost functions
# Tobler1993 velocity
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

## ## plot cost function
plot(tobler1993a
     ,xlim = c(-1,1)
     ,main = "Slope-dependent cost function (Tobler)"
     ,xaxt="n",yaxt="n"
     ,xlab = "",ylab = "")
mtext(side = 1, text = "slope", line = 2)
mtext(side = 2, text = "speed (km/h)", line = 2)
axis(2, mgp=c(3, .5, 0))
axis(1, mgp=c(3, .5, 0))
abline(v = 0, lty = 2)

library(gdistance) # https://cran.r-project.org/web/packages/gdistance/index.html
# auxilliary function
hdiff <- function(x){(x[2]-x[1])} # calculates the difference of a vector

# transitional object 
hd <- transition(ras,hdiff,nz,symm=FALSE) # calculate the slope of the input data; store the conductivity for the movement to all other cells from any specfic cell
## Warning in .TfromR(x, transitionFunction, directions, symm): transition
## function gives negative values
slope <- geoCorrection(hd,scl=FALSE) # geocorrection of the object because the diagonal distance is larger than a vertical distance [in the pixel case]
adj <- adjacent(x=ras, cells=1:ncell(ras), direction=nz) # adjacency list - from X to Y
cost <- slope       
cost[adj] <- tobler1993a(slope[adj]) # and here the conductivity function comes into account
conduct <- geoCorrection(cost, scl=FALSE) # transform cost to conductivity; conductivity=cost/dist; time=1/conductivity; we need the geocorection twice because


p.colombo <- spTransform(p.colombo, CRSobj = srtm@crs)
p.kandy <- spTransform(p.kandy, CRSobj = srtm@crs)

co.to.ka <- shortestPath(conduct, origin = p.colombo@coords, goal = p.kandy@coords, output="SpatialLines")
ka.to.co <- shortestPath(conduct, origin = p.kandy@coords, goal = p.colombo@coords, output="SpatialLines")


plot(srtm,
     xlim = c(300000,500000), ylim = c(730000,830000),
     main = "Least-cost path\n based on Tobler's function for walking speed")
lines(co.to.ka, col = "red")
lines(ka.to.co, col = "blue")
points(p.colombo, pch = 19, cex = .5)
points(p.kandy, pch = 19, cex = .5)
text(p.colombo@coords[1], p.colombo@coords[2] - 2500, "Colombo")
text(p.kandy@coords[1], p.kandy@coords[2] + 2500, "Kandy")

## insert the modern highway from Colombo to Kandy as comparison
m.r <- readOGR("./data/DSL250-Shp", "Roads")
## OGR data source with driver: ESRI Shapefile 
## Source: "./data/DSL250-Shp", layer: "Roads"
## with 3358 features
## It has 10 fields
m.r <- spTransform(m.r, srtm@crs)
lines(m.r[m.r@data$CLASS=="A1",], col="black")

## and at the the a nice legend
legend("bottom", c("Colombo - Kandy", "Kandy - Colombo", "A1"),
       lty = c(1,1,1), col = c("red", "blue", "black"), bty = "n")

e.co.ka <- SpatialLines(list(Lines(Line(coords = (rbind(p.colombo@coords,p.kandy@coords))), ID="1")),proj4string = srtm@crs)

p.co.to.ka <- extract(x=srtm, y=co.to.ka, along=TRUE)
p.ka.to.co <- extract(x=srtm, y=ka.to.co, along=TRUE)
p.e.co.ka <- extract(x=srtm, y=e.co.ka, along=TRUE)

library(rgeos)

par(mfrow = c(3,1))
plot(p.e.co.ka[[1]], type="l",
     main = "Elevation profile \nfrom Colombo to Kandy (euclidean)",
     sub = paste("Distance ",round(gLength(e.co.ka)/1000,2)," km",sep=""))
plot(p.co.to.ka[[1]], type="l",
     main = "Elevation profile \nfrom Colombo to Kandy (least cost)",
     sub = paste("Distance ",round(gLength(co.to.ka)/1000,2)," km",sep=""))
plot(p.ka.to.co[[1]], type="l",
     main = "Elevation profile \nfrom Kandy to Colombo (least cost)",
     sub = paste("Distance ",round(gLength(ka.to.co)/1000,2)," km",sep=""))