# (c) 2016 Dr. Karsten Schmidt, M. Sc. Philipp M. Gries
# Collaborative Research Centre 1070 - ResourceCultures
# eScience Centre - University of Tübingen


############
# Packages #
############

options(repos='http://cran.rstudio.com/')
install.packages("raster")
install.packages("gdistance")
install.packages("sp")
install.packages("lattice")
install.packages("gstat")
install.packages("rgeos")

######################
# PART 1 - Libraries #
######################

library(raster)
library(gdistance)
library(sp)
library(lattice)
library(gstat)
library(rgeos)

#####################################
# PART 2 - Tobler's Hiking Function #
#####################################

# implementing tobler's hiking function
ToblersHikingFunction <- function(x){ 6 * exp(-3.5 * abs(tan(x*pi/180) + 0.05)) }

# plot tobler's hiking function
TheoreticalSlopes <- seq(-70,70,1)
WlkSpeed       <- ToblersHikingFunction(TheoreticalSlopes)
plot(TheoreticalSlopes, WlkSpeed, type="l", col ="red", lwd = 2, lty="dashed",  ylab="walking speed [km/hr]", xlab="Average slope in degrees", axes=F)
axis(1, tck=-.01, at= TheoreticalSlopes[seq(1,length(TheoreticalSlopes),10)], labels= TheoreticalSlopes[seq(1,length(TheoreticalSlopes),10)])
axis(2)
abline(v=0, lty="dashed", col ="gray")
title(expression("The Tobler Hiking Function\nspeed = 6 * exp(-3.5 * abs(Slope + 0.05)"))

####################
# PART 3 - Setting #
####################

# 1 DIRECTORIES
inDir <- "C:/PATH/TO/YOUR/DATA"  # directory to read data
outDir <- "C:/PATH/TO/YOUR/DATA" # directory to write data

# 2 DATASETS
# filenames of input data and x-y-location of initial point
DEM = "SET_SRTM90.asc"
SLOPE = ""
POINT <- c(3437932,5332185) # c(X-coordinate,Y-coordinate), e.g. c(1234567,1234567)

# enter file names and datatype of output data 
TCR = "TimeCostRasterExample0802" # name of time cost raster
SLG = "SlopeGradientExample0802" # name of slope gradient raster
CTL = "ContourLinesExample0802" # name of contour lines shapefile
rdt = "ascii" # datatype of output raster files, for more information use 'help(writeRaster)'

# 3 ENVIRONMENTAL VARIABLES
# slope
# only if no SLOPE input raster available
# select algorithm to calculate slope
# 4 == 4 neighbours, smooth surfaces, Fleming and Hoffer (1979), Ritter (1987) (see raster package terrain)
# 8 == 8 neighbours, rough surfaces, Horn (1981) (see raster package terrain)
numberOfNeighbors <- 8

# damp hiking speed if slope is steeper than the dampingFactor (slope in degrees)
damping <- FALSE
dampingFactor = 16 

# move cases for transition layer
# select move case for calculating time-cost
# cell connections, directions for transition:
# 4 == 4 directions, Rook case
# 8 == 8 directions, Queen case
# 16 == 16 directions, knight move
numberOfDirections <- 8

# time interval of interest
# select a time interval of interest:
timeOfInterest <- 2 # use integer > 0
#timeOfInterest <- -1 #  for no specific time interval
# Warning: If -1, time-cost-analysis of complete input raster is used (full extent) 
# and result in sesquipedalian computing time and/or may cause memory errors

# select a number and interval of isochrones output:
numberOfIsochrones <- 2 # number of isochrones [hours], integer
intervallOfIsochrones <- 1 # intervals of insochrones [hours], numeric

#####################################################
# PART 4 - Reading data and create spatial datasets #
#####################################################

setwd(inDir)

# ELEVATION
rDEM <- raster(DEM) # convert DEM to raster datatype

# SLOPE GRADIENT
rSLOPE <- NULL # empty slope raster variable
if (nchar(SLOPE) > 0) { # if SLOPE is set
  rSLOPE <- raster(SLOPE) # convert SLOPE to raster datatype
} else { # if SLOPE is not set
  # calculate SLOPE from DEM
  if(!is.na(projection(rDEM))) {
    rSLOPE <- terrain(rDEM, opt='slope', unit='degrees', neighbors=numberOfNeighbors) # calculate slope from elevation model
  } else {
    # no projection of elevation model result in projection error. Create projection file. 
    print("PROJECTION ERROR: no projection is set for ELEVATION input file.")
  }
}

# create SPATIAL POINT
SPATIALPOINT <- data.frame(x=POINT[1],y=POINT[2])
coordinates(SPATIALPOINT) <- ~ x+y
projection(SPATIALPOINT) <- projection(rDEM)

# clip raster if time limit is set
rSLOPE4TimeCost <- rSLOPE
rDEM4Statistics <- rDEM

if (timeOfInterest > 0) {
  
  # maximal distance in limit + time buffer 0.25h
  maxHikingDistance <- round(max(WlkSpeed)*(timeOfInterest+0.25)*1000)
  
  # buffer hiking distance around point of interest
  bufferMaxHikingDistance <- buffer(SPATIALPOINT, maxHikingDistance)
  
  #clip raster and set new extent
  rDEM_clip <- crop(rDEM,extent(bufferMaxHikingDistance))
  rSLOPE_clip <- crop(rSLOPE,extent(bufferMaxHikingDistance))
  
  # DEM raster for statistic section below
  rDEM4Statistics <- rDEM_clip
  
  # slipped slope raster for time cost calculation
  rSLOPE4TimeCost <- rSLOPE_clip
}

################################
# PART 5 -  Time Cost Analysis #
################################

# Calculate velocity per cell (V)
rVelocity.kmh <- calc(rSLOPE4TimeCost, ToblersHikingFunction) # s wird ab hier durch slope4TimeCost ersetzt

# Velocity in m/s
rVelocity.ms <- calc(rVelocity.kmh, fun=function(x) { ((x*1000)/3600) })

# damping velocity
# damp the hiking speed if slope is greater than the dampingFactor
# default: from a slope grater than 16 degrees the hiking speed is reduced by extremly exaggerating the slope 
# if damping is true
if (damping) {
  rDamping <- rSLOPE4TimeCost
  rDamping[rDamping >  dampingFactor]  = 1000
  rDamping[rDamping <= dampingFactor]  = 1  
  rVelocity.ms <- rVelocity.ms/rDamping

}

# transition layers
# Rook move [4]; Queens move [8]; Knights move [16]
lTransition <- transition(rVelocity.ms, transitionFunction=mean, directions=numberOfDirections)

# With GeoCorrection
# take diagonal and cardinal connections into account
lGeoCorrection <- geoCorrection(lTransition, type="r")
rAccumulatedCostSurface.s <- accCost(lGeoCorrection, SPATIALPOINT)

# convert time [sec] to [hour]
rAccumulatedCostSurface.h <- rAccumulatedCostSurface.s/3600

# plot results
plot(rAccumulatedCostSurface.h)
plot(SPATIALPOINT,add=TRUE)

# create contour lines
vContourLines <- rasterToContour(rAccumulatedCostSurface.h, levels=seq(0, numberOfIsochrones, intervallOfIsochrones))
#vContourLines <- rasterToContour(rAccumulatedCostSurface.h, levels=seq(0, 2, 0.25))
plot(vContourLines, add=TRUE)

#############################
# PART 6 - Zonal Statistics #
#############################

# create dataframes for statistics
zonalStatistics   <- data.frame(matrix(0,2,5))
statNames <- c('1st hour','min','max','mean','sd')
names(zonalStatistics) <- statNames
zonalStatistics[1,1] <- 'DEM'
zonalStatistics[2,1] <- 'SLOPE'

# create zones
rasterZones <- rAccumulatedCostSurface.h
rasterZones[rAccumulatedCostSurface.h <= 1] = 1 # pixels of first hour
rasterZones[rAccumulatedCostSurface.h > 1]  = 2 # pixels more than first hour

# DEM, SLOPE statistics first hour
for(st in 2:length(statNames)) {
    zonalStatistics[1,st] <- zonal(rDEM4Statistics, rasterZones, statNames[st])[1,2]
    zonalStatistics[2,st] <- zonal(rSLOPE4TimeCost, rasterZones, statNames[st])[1,2]
} 

print(zonalStatistics)
                    
########################
# PART 7 - Write Files #
########################

setwd(outDir) # set directory of output files
writeRaster(rAccumulatedCostSurface.h,TCR, format=rdt) # write time cost raster
if (nchar(SLOPE) == 0) {
  writeRaster(rSLOPE4TimeCost,SLG, format=rdt) # write slope raster
}
shapefile(vContourLines,filename=CTL) # write contourlines shapefile