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)
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
## 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=""))