## ----setup, include=FALSE-----------------------------------------------------
setwd("C:/github/reagro/web/source/cases/_R")
library(gdistance)


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


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


## ----gadm, message=FALSE, warning=FALSE---------------------------------------
adm <- gadm("TZA", level = 1, path=".")
roads <- geodata::osm("Tanzania", "highways", ".")


## ----simplr, eval=FALSE-------------------------------------------------------
## roads <- simplifyGeom(roads, tolerance=0.01)


## ----plotRoads----------------------------------------------------------------
plot(tza_slope)
lines(roads)
lines(roads[roads$highway == "secondary", ], lwd=2, col="blue")
lines(roads[roads$highway == "primary", ], lwd=4, col="red")


## ----rasterize----------------------------------------------------------------
cfile <- "rdcost.tif"
roadtypes <- c("primary", "secondary", "tertiary")

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, field=roads$speed, filename=cfile, wopt=list(names="slope_cost"), overwrite=TRUE)

	} else {
	# no need to repeat if we already have done this
	rd_cost <- rast(cfile)
} 

# aggregate to improve visualization
a <- aggregate(rd_cost, 3, min, na.rm=TRUE)
plot(a, col=c("black", "blue", "red"), main="Travel cost (min/m)")


## -----------------------------------------------------------------------------
tza_lc <- agrodata::reagro_data("TZA_globcover")
plot(tza_lc, main = "landcover classes")
lines(adm)


## ---- echo = FALSE------------------------------------------------------------
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"),
  Travel_speed=c(0.04, 0.04, 0.04, 0.03, 0.05, 0.01, 0.01, 0.11, 0.13)
)

knitr::kable(text_tbl)


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

#reclassifying
tza_lc_cost <- classify(tza_lc, rc) 


## -----------------------------------------------------------------------------
lcfname <- "lc_cost.tif"
if (!file.exists(lcfname)) {
  # first aggregate to about the same spatial resolution
  lc_cost <- aggregate(tza_lc_cost, 3, mean)
  # then resample
  lc_cost <- resample(lc_cost, tza_slope, filename=lcfname, wopt=list(names="lc_cost"), overwrite=TRUE)
} else {
  lc_cost <- rast(lcfname)
}


## -----------------------------------------------------------------------------
plot(lc_cost, main = "Off-road travel costs (min/m) based on land cover class")


## -----------------------------------------------------------------------------
# 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)
cost <- cost * slope_cost

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


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


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

cost <- raster(cost)
conductance <- 1/cost

#Creating a transition object 
tran <- transition(conductance, transitionFunction=mean, directions= 8)


## -----------------------------------------------------------------------------
tran <- geoCorrection(tran, type="c")


## -----------------------------------------------------------------------------
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 <- cbind(lon, lat)
spCities <- sp::SpatialPoints(cities)

#Estimating 
Ac <- accCost(tran, fromCoords=spCities)


## ----costrast-----------------------------------------------------------------
A <- rast(Ac) / 60
AA <- clamp(A, 0, 24) |> mask(tza_slope)
plot(AA, main="Access to markets in Tanzania ()")
lines(roads)
points(cities, col="blue", pch=20, cex=1.5)