## ----setup, include=FALSE-----------------------------------------------------


## ----slope, message=FALSE-----------------------------------------------------
library(geodata)
tza_alt <- elevation_30s("Tanzania", path=".")
tza_slope <- slope(tza_alt, unit="radians", neighbors=8)
plot(tza_slope, main="Slope")


## ----access20-----------------------------------------------------------------
decay <- 1.5
slope_cost <- exp( decay * tan(tza_slope) )
names(slope_cost) <- "slope_cost"


## ----access30, message=FALSE, warning=FALSE-----------------------------------
adm <- gadm("TZA", level = 1, path=".")


## ----access40-----------------------------------------------------------------
library(osmdata)

getOSMlines <- function(obj, key, value) {
    b <- sp::bbox(obj)  #returns the boundary box of the object
    q <- osmdata::opq(b)  #Creates an Overpass query within the boundary box
    q <- osmdata::add_osm_feature(q, key, value)  #queries the key and value atrributes
    d <- osmdata::osmdata_sp(q)  #converts the data to a sp object
    d$osm_lines[, key]   #returns only the lines found in the search
}
roadtypes <- c("primary", "secondary", "tertiary")



## ----access_getOSM, eval=FALSE------------------------------------------------
## x <- lapply(seq_along(adm), function(i) getOSMlines(adm[i, ], "highway", roadtypes))


## ----access_processOSM, eval=FALSE--------------------------------------------
## # combine the results for each area
## rd <- do.call(bind, x)
## # remove anything outside of Tanzania
## rd <- crop(rd, adm)
## # aggregate by type, also removing duplicates line parts
## roads <- aggregate(rd, "highway")


## ----access_simplr, eval=FALSE------------------------------------------------
## library(rmapshaper)
## roads <- ms_simplify(roads, keep=0.01)


## ----access_plotRoads---------------------------------------------------------
library(reagrodata)
roads <- reagro_data("tza_roads")
plot(tza_slope)
lines(roads)
lines(roads[roads$highway == "secondary", ], lwd=2, col="blue")
lines(roads[roads$highway == "primary", ], lwd=4, col="red")


## ----access_rasterize---------------------------------------------------------
cfile <- "rdcost.tif"

if (!file.exists(cfile)) {

	i <- match(roads$highway, roadtypes)
	roads$speed <- c(0.001, 0.0015, 0.002)[i]
	
	rd_cost <- rasterize(roads, tza_slope, values="speed", filename=cfile, wopt=list(names="slope_cost"))
	
} else {
	# no need to repeat if we already have done this
	rd_cost <- rast(cfile)
} 

plot(rd_cost, main="Travel cost (min/m")


## ----access_lc----------------------------------------------------------------
tza_lc <- reagro_data("TZA_globcover")
plot(tza_lc, main = "GLOBCOVER 2009 v2.3 land classes")
lines(adm)


## ----accesstab, echo = FALSE--------------------------------------------------
library(knitr)
text_tbl <- data.frame(
  Value = c(40,50,70,160,170,190,200,210,220),
  LandClass = c("Closed to open (>15%) broadleaved evergreen or semi-deciduous forest (>5m)",
                "Closed (>40%) broadleaved deciduous forest (>5m)",
                "Closed (>40%) needleleaved evergreen forest (>5m)",
                "Closed to open (>15%) broadleaved forest regularly flooded (semi-permanently or temporarily) - Fresh or brackish water",
                "Closed (>40%) broadleaved forest or shrubland permanently flooded - Saline or brackish water",
                "Artificial surfaces and associated areas (Urban areas >50%)",
                "Bare areas",
                "Water bodies",
                "Permanent snow and ice"),
  TravelCost=c(0.04, 0.04, 0.04, 0.03, 0.05, 0.01, 0.01, 0.11, 0.13)
)

kable(text_tbl)


## ----access200----------------------------------------------------------------
rc <- data.frame(from=as.vector(unique(tza_lc)), to=0.02)
rc$from[rc$from %in% c(190,200)] <- 0.01
rc$from[rc$from == 160] <- 0.03
rc$from[rc$from %in% c(40,50,70)] <- 0.04
rc$from[rc$from == 170] <- 0.05
rc$from[rc$from == 210] <- 0.11
rc$from[rc$from == 220] <- 0.13

#reclassifying
tza_lc_cost <- classify(tza_lc, rc) 


## ----access210----------------------------------------------------------------
lcfname <- "lc_cost.tif"
if (!file.exists(lcfname)) {
  lc_cost <- warp(tza_lc_cost, tza_slope, filename=lcfname, wopt=list(names="lc_cost"))
} else {
  lc_cost <- rast(lcfname)
}


## ----access220----------------------------------------------------------------
plot(lc_cost, main = "Travel costs (min/m) depending on land cover class")


## ----access230----------------------------------------------------------------
# Combine the cost layers
all_cost <- c(rd_cost, lc_cost)

#getting the minimum value of each grid cell
cost <- min(all_cost, na.rm=TRUE)*slope_cost

writeRaster(cost, filename = "cost.tif", overwrite=TRUE)

plot(cost, main="Final cost layer (min/m)")


## ----access240, eval=FALSE----------------------------------------------------
## install.packages("gdistance")


## ----access250----------------------------------------------------------------
# Combine the cost layers
library(gdistance)

cost <- raster("cost.tif")

#Creating a transition object 
trans <- transition(cost,transitionFunction=mean, directions= 8)




## ----access260----------------------------------------------------------------
trans <- geoCorrection(trans, type="c")


## ----access270----------------------------------------------------------------
lat=c(-6.17, -6.81, -5.02, -2.51, -3.65, -8.90, -3.34, -3.36, -10.67)
lon=c(35.74, 37.66, 32.80, 32.90, 33.42, 33.46, 37.34, 36.68, 35.64)
cities <- SpatialPoints(cbind(lon, lat))

#Estimating 
A <- accCost(trans, fromCoords=cities)

plot(A, main="Access to markets in Tanzania ()")
points(cities)