# empty memory and workspace
gc()
rm(list=ls())

# check directory
getwd()

### to install all required R packages
ls <- c("sp", "raster", "rgdal", "mapview", "corrplot",
        "caret", "rpart", "rpart.plot", "randomForest")
new.packages <- ls[!(ls %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

# load library
library(sp)
library(raster)
library(rgdal)
library(mapview)
library(corrplot)
library(caret)
library(rpart)
library(rpart.plot)
library(randomForest)

# import raster layers of covariates (Remote Sensing data)
blue = raster("./RS/Blue.tif")
green = raster("./RS/Green.tif")
red = raster("./RS/Red.tif")
nir = raster("./RS/Nir.tif")

# calculate some Remote Sensing indices such as NDVI
ndvi = (nir - red) / (nir + red)

# plot raster layers of covariates (Remote Sensing data)
plot(blue, main = "blue band of landsat", xlab = "Easting (m)", ylab = "Northing (m)")

plot(red, main = "red band of landsat", xlab = "Easting (m)", ylab = "Northing (m)")

plot(ndvi, main = "NDVI of landsat", xlab = "Easting (m)", ylab = "Northing (m)")

# make a stack layer from Remote Sensing data
landsat = stack(blue,green,red,nir,ndvi)

# inspect the structure of stack layer
landsat

# set the names to the rasters in the stack file
names(landsat) <- c("blue","green","red","nir","ndvi")

# plot the stack layers of Remote Sensing data
plot(landsat)

# import the terrain data using list.files command
dem_lst <- list.files("./Terrain/", pattern="\\.sdat$", full.names = TRUE)

# inspect the structure of list
dem_lst

# make a stack layer from terrain data 
terrain = stack(dem_lst)

# plot the stack layers of terrain data
plot(terrain)

# Re-sampling Remote Sensing data based on terrain data
RS_terrain = resample(landsat,terrain, method="bilinear")

# make a stack file from re-sampled remote sensing data and terrain data
covariates_stack = stack(terrain, RS_terrain)

# inspect the structure of stack layer
names(covariates_stack)

# import shape file of point data
point = readOGR("./Shapefile","points")

# import shape file of Bavaria boundary
Bavaria = readOGR("./Shapefile","Bavaria")

# plots of the shape files
plot(Bavaria)
plot(point, add=T)

# plot an interactive map
mapview::mapview(point, zcol = "OC", at = c(0,5,10,15,20,25,30,35,200), legend = TRUE)

# plot the point on the raster
plot(covariates_stack$DEM, main = "DEM", xlab = "Easting (m)", ylab = "Northing (m)")
plot(point, add =T,pch = 19)
plot(Bavaria, add =T)

# inspect the point shape files
# types of the spatial data
class(point)

# summary statistics of data
summary(point)

# histogram of data
hist(point$OC,col ="blue", xlab= "SOC (g/kg)", main="Histogram")

# extract covariate values at each point of observation 
cov = extract(covariates_stack,point,method='bilinear', df=TRUE)

# inspect the data.frame of cov
str(cov)

# combining covariates and soil properties
cov_soil = cbind(cov[,-1], OC=point$OC)

# inspect the data.frame of cov_soil
str(cov_soil)

# check the correlation covariates and OC
corrplot.mixed(cor(cov_soil), lower.col = "black", number.cex = .7)

# remove na values
cov_soil <- cov_soil[complete.cases(cov_soil),]

# remove high values of oc
cov_soil <- cov_soil[cov_soil$OC<78,]

# check number of rows
nrow(cov_soil)

# check number of column 
ncol(cov_soil)

# split the data to training (70%) and testing (30%) sets
trainIndex <- createDataPartition(cov_soil$OC, 
                                  p = 0.7, list = FALSE, times = 1)

# subset the datasets
cov_soil_Train <- cov_soil[ trainIndex,]
cov_soil_Test  <- cov_soil[-trainIndex,]

# inspect the two datasets
str(cov_soil_Train)

str(cov_soil_Test)  

# fit models on training set: regression model
linear_fit <- lm(OC ~ CNBL+DEM+MCA+Slope+TWI+blue+green+red+nir+ndvi, 
                 data=cov_soil_Train)

# look at the summary of linear model
summary(linear_fit)

# apply the linear model on testing data
OC_linear_Pred <- predict(linear_fit, cov_soil_Test)  

# check the plot actual and predicted OC values
plot(cov_soil_Test$OC, OC_linear_Pred, main="Linear model", 
     col="blue",xlab="Actual OC", ylab="Predicted OC", 
     xlim=c(0,100),ylim=c(0,100))
abline(coef = c(0,1),  col="red" )

# calculate correlation
cor_linear <- cor(cov_soil_Test$OC, OC_linear_Pred)
cor_linear

# calculate RMSE
RMSE_linear <- sqrt(mean((cov_soil_Test$OC - OC_linear_Pred)^2))
RMSE_linear

# fit models on training set: decision tree model
tree_fit <- rpart(OC ~ CNBL+DEM+MCA+Slope+TWI+blue+green+red+nir+ndvi,
                  method="anova", data=cov_soil_Train,
                  control = rpart.control(minsplit = 15, cp = 0.0001))

# display the results
printcp(tree_fit)

# visualize cross-validation results
plotcp(tree_fit) 

# detailed summary of splits
summary(tree_fit) 

# visualize the tree
rpart.plot(tree_fit) 

# apply the tree model on testing data
OC_tree_Pred <- predict(tree_fit, cov_soil_Test)

# check the plot actual and predicted OC values
plot(cov_soil_Test$OC, OC_tree_Pred, main="Tree model", 
     col="blue",xlab="Actual OC", ylab="Predicted OC", xlim=c(0,100),ylim=c(0,100))
abline(coef = c(0,1),  col="red" )

# calculate correlation
cor_tree <- cor(cov_soil_Test$OC, OC_tree_Pred)
cor_tree

# calculate RMSE
RMSE_tree <- sqrt(mean((cov_soil_Test$OC - OC_tree_Pred)^2))
RMSE_tree

# fit models on training set: random forest model
rf_fit <- randomForest(OC ~ CNBL+DEM+MCA+Slope+TWI+blue+green+red+nir+ndvi, 
                       data=cov_soil_Train,ntree=1000, do.trace = 25)

# detailed summary of random forest
summary(rf_fit) 

# visualize the importance of random forest
varImpPlot(rf_fit) 

# apply the random forest model on testing data
OC_rf_Pred <- predict(rf_fit, cov_soil_Test)

# check the plot actual and predicted OC values
plot(cov_soil_Test$OC, OC_rf_Pred, main="Tree model", 
     col="blue",xlab="Actual OC", ylab="Predicted OC", xlim=c(0,100),ylim=c(0,100))
abline(coef = c(0,1),  col="red" )

# calculate correlation
cor_rf <- cor(cov_soil_Test$OC, OC_rf_Pred)
cor_rf

# calculate RMSE
RMSE_rf <- sqrt(mean((cov_soil_Test$OC - OC_rf_Pred)^2))
RMSE_rf

# check the accuracy of three models
RMSE_models <- c(Linear=RMSE_linear,Tree=RMSE_tree,RF=RMSE_rf)
cor_models <- c(Linear=cor_linear,Tree=cor_tree,RF=cor_rf)

# plot the final results
par(mfrow=c(1,2))
barplot(RMSE_models, main="RMSE",col=c("red","blue","green"))
barplot(cor_models, main="Correlation",col=c("red","blue","green"))

# apply the best model on the stack layer
map_rf <- raster::predict( covariates_stack,rf_fit)

# plot the final maps
spplot(map_rf, main="SOC map based on RF model")

# plot an interactive map
mapview::mapview(map_rf)