library(leaflet)
library(sf)
# Call external functions
source(paste0( scripts_path , "VMCPTF_function.R")) #
source(paste0( scripts_path , "EffPAWHC_function.R")) #
source(paste0( scripts_path , "RI_function.R")) #
source(paste0( scripts_path , "EvaluateRI_function.R")) #
source(paste0( scripts_path , "Spline_function.R")) #
source(paste0( scripts_path , "Auxiliary_functions.R")) #
## ----read_aoi, collapse=T, results='hide', message=F--------------------------------------------------------
# read a shapefile with the area of interest
aoi_read <- sf::st_read(paste0("//WURNET.NL/DFS-Root/PROJECTS/ISRIC_Workspace/projects/5318018099_S2P/DST_prototype/model/layersMisc/admin/admin1_madagascar_gadm.shp"))
aoi_select <- aoi_read[aoi_read$NAME_1 == "Antananarivo", ]
aoi <-  st_geometry(aoi_select)
# List maps names
maps_lst <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids")
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
# Exchangeable sodium (Na+) "ExchNa" layers available in 1km
sd_lst_2 = c("_sd1_1km.tif","_sd2_1km.tif","_sd3_1km.tif","_sd4_1km.tif","_sd5_1km.tif","_sd6_1km.tif")
for (p in 1:length(sd_lst_2)) {
exna = "gyga_af_enax_t__m"
print(paste0(exna, sd_lst_2[p]))
name_tmp = rast(paste0(working_path,exna,sd_lst_2[p]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign("ExchNa",stack_tmp)
names("ExchNa")
View(ExchNa)
ExchNa
names(ExchNa)
View(ExchAcids)
for (p in 1:length(sd_lst_2)) {
raster_tmp = list()
exna = "gyga_af_enax_t__m"
print(paste0(exna, sd_lst_2[p]))
name_tmp = rast(paste0(working_path,exna,sd_lst_2[p]))
names(name_tmp) <- paste0("ExchNa",sd_lst[p])
update(name_tmp, names=TRUE)
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign("ExchNa",stack_tmp)
names(ExchNa)
ExchNa
p = 1
raster_tmp = list()
exna = "gyga_af_enax_t__m"
print(paste0(exna, sd_lst_2[p]))
name_tmp = rast(paste0(working_path,exna,sd_lst_2[p]))
name_tmp
names(name_tmp) <- paste0("ExchNa",sd_lst[p])
names(name_tmp
)
update(name_tmp, names=TRUE)
name_tmp
raster_tmp = c(raster_tmp,name_tmp)
raster_tmp
p = 2
raster_tmp = list()
for (p in 1:length(sd_lst_2)) {
exna = "gyga_af_enax_t__m"
print(paste0(exna, sd_lst_2[p]))
name_tmp = rast(paste0(working_path,exna,sd_lst_2[p]))
names(name_tmp) <- paste0("ExchNa",sd_lst[p])
update(name_tmp, names=TRUE)
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign("ExchNa",stack_tmp)
ExchNa
names(ExchNa)
Sand
Sand[[2]]
UpDpth <-c(0,5,15,30,60,100)
LowDpth <-c(5,15,30,60,100,200)
for (i in 1:length(maps_lst)) {
print (kable(stats.all(get(maps_lst[i]))))
print(maps_lst[i])
}
maps_lst2 <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids",
"ExchNa")
for (i in 1:length(maps_lst2)) {
print (kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
}
maps_lst2 <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids",
"ExchNa")
for (i in 1:length(maps_lst2)) {
print (kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
}
length(maps_lst2)
maps_lst <- c(#"Sand",
#"Silt",
#"Clay",
#"CrsFrg",
#"OrgC",
"BulkDens"#,
#"CEC" ,
#"pHH2O",
#"ExchBases",
#"EC",
#"ExchAcids"
)
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
print("Done")
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
print("Done")
"a"
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
print("Done")
i = 5
print (kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
i = 6
print(kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
BulkDens
maps_lst <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
#"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids")
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
BulkDens
OrgC
CEC
maps_lst <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
#"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids")
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
names(name_tmp) <- paste0(maps_lst[i],sd_lst[p])
update(name_tmp, names=TRUE)
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
BulkDens
i = 6
print(kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
maps_lst2 <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids",
"ExchNa")
for (i in 1:length(maps_lst2)) {
print(kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
}
Sand
i = 1
p = 1
raster_tmp = list()
j = 1
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
name_tmp
paste0(working_path,maps_lst[i], sd_lst[j])
raster_tmp = list()
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
name_tmp
paste0(maps_lst[i],sd_lst[p])
names(name_tmp) <- paste0(maps_lst[i],sd_lst[p])
update(name_tmp, names=TRUE)
update(name_tmp, names=TRUE)
raster_tmp = c(raster_tmp,name_tmp)
raster_tmp
paste0(maps_lst[i])
maps_lst <- c("Sand",
"Silt",
"Clay",
"CrsFrg",
"OrgC",
"BulkDens",
"CEC" ,
"pHH2O",
"ExchBases",
"EC",
"ExchAcids")
# List sd's
sd_lst = c("_sd1.tif","_sd2.tif","_sd3.tif","_sd4.tif","_sd5.tif","_sd6.tif")
# Each map has 6 sd's
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
names(name_tmp) <- paste0(maps_lst[i],sd_lst[p])
update(name_tmp, names=TRUE)
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
Sand
maps_lst[i]
sd_lst[p]
for (i in 1:length(maps_lst)) { # all the input raster layers are read using this loop
raster_tmp = list()
for (j in 1:length(sd_lst)) {
print(paste0(maps_lst[i], sd_lst[j]))
name_tmp = rast(paste0(working_path,maps_lst[i], sd_lst[j]))
names(name_tmp) <- paste0(maps_lst[i],sd_lst[j])
update(name_tmp, names=TRUE)
raster_tmp = c(raster_tmp,name_tmp)
}
stack_tmp = rast(raster_tmp)
assign(paste0(maps_lst[i]),stack_tmp)
}
for (i in 1:length(maps_lst2)) {
print(kable(stats.all(get(maps_lst2[i]))))
print(maps_lst2[i])
}
all.properties.sd <- read.data(paste0(maps_lst2, "_sd1.tif"))
all.properties.sd <- paste0(maps_lst2, "_sd1")
plot.map(all.properties.sd, fixed=FALSE)
all.properties.sd <- read.data(paste0(maps_lst2, "_sd2.tif"))
all.properties.sd2 <- paste0(maps_lst2, "_sd1")
prueba = ExchNa["ExchNa_sd5"]
t = 1
paste0(maps_lst2[t],"_sd1")
band_tmp = maps_lst2[t][paste0(maps_lst2[t],"_sd1")]
maps_lst2[t]
band_tmp = maps_lst2[t]
band_tmp = Sand
band_tmp = list(Sand,
Silt,
Clay,
CrsFrg,
OrgC,
BulkDens,
CEC ,
pHH2O,
ExchBases,
EC,
ExchAcids,
ExchNa)
rast_lyrs = list(Sand,
Silt,
Clay,
CrsFrg,
OrgC,
BulkDens,
CEC ,
pHH2O,
ExchBases,
EC,
ExchAcids,
ExchNa)
rast_lyrs = list(Sand,
Silt,
Clay,
CrsFrg,
OrgC,
BulkDens,
CEC ,
pHH2O,
ExchBases,
EC,
ExchAcids,
ExchNa)
t = 1
band_tmp = rast_lyrs[paste0(maps_lst2[t],"_sd1")]
paste0(maps_lst2[t],"_sd1")
band_tmp = rast_lyrs[t]
band_tmp
prueba = rast_lyrs[1]["Sand_sd1"]
rast_lyrs[1]["Sand_sd1"]
rast_lyrs[1]
rast_lyrs[[1]]
rast_lyrs[[1]]["Sand_sd1.tif"]
rast_lyrs[[1]][1]
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
band_tmp
for (t in 1:length(sd_lst)){
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
}
for (t in 1:length(sd_lst)){
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
raster_tmp = c(raster_tmp,band_tmp)
}
raster_tmp = list()
for (t in 1:length(sd_lst)){
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
raster_tmp = c(raster_tmp,band_tmp)
}
raster_tmp = list()
for (t in 1:length(sd_lst)){
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
raster_tmp = c(raster_tmp,band_tmp)
} stack_sd1 = rast(raster_tmp)
raster_tmp = list()
for (t in 1:length(sd_lst)){
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
raster_tmp = c(raster_tmp,band_tmp)
}
stack_sd1 = rast(raster_tmp)
plot.map(stack_sd1, fixed=FALSE)
stack_sd1
raster_tmp = list()
for (t in 1:length(rast_lyrs)){
band_tmp = rast_lyrs[[t]][paste0(maps_lst2[t],"_sd1")]
raster_tmp = c(raster_tmp,band_tmp)
}
stack_sd1 = rast(raster_tmp)
stack_sd1
plot.map(stack_sd1, fixed=FALSE)
paste0(working_path,  "agland_crop.tif")
mask=rast(paste0(working_path,  "agland_crop.tif"))
?app
mask = app(mask, function(x) {ifelse(is.na(x),  NA, x)} )
mask
pHH2O[[1]]
mask = app(mask, function(x) {ifelse(is.na(x),  NA, 1)} )
plot(mask)
drainClasses=rast(paste0(working_path,  "drainClasses.tif")) # read the drain classes map
pHH2O <- pHH2O /10
sum.tex <- Clay+Silt+Sand
Clay <- Clay/(sum.tex)*100
Silt <- Silt/(sum.tex)*100
Sand <- Sand/(sum.tex)*100
sum.tex <- Clay+Silt+Sand
VMCPTF <- list(NULL)
for (x in 1:6) { # loop to apply the function to 6 depths
# compile in a raster stack all the input data needed for the same layer depth
input <- c(Sand[[x]] , Silt[[x]], Clay[[x]], OrgC[[x]], BulkDens[[x]], CEC[[x]], pHH2O[[x]], CrsFrg[[x]],ExchBases[[x]], EC[[x]], ExchAcids[[x]])
# adapt the names of the raster compiled layers
names(input) <- c("Sand", "Silt", "Clay",   "OrgC",  "BulkDens", "CEC" , "pHH2O", "CrsFrg", "ExchBases" , "EC", "ExchAcids")
# apply the function VMCPTF.function
VMCPTF[[x]] <- VMCPTF.function( Sand=input[["Sand"]],   Silt=input[["Silt"]], Clay=input[["Clay"]],
OrgC=input[["OrgC"]], BulkDens=input[["BulkDens"]], CEC=input[["CEC" ]],
pHH2O=input[["pHH2O"]],  h1=-10, h2=-20, h3=-31.6, h4=-1585,
PTF.coef, fix.values=TRUE, print.coef=TRUE)
} # h0 is defined and then VMC_h0 is calculated in the function
# compile in a raster stack all the input data needed for the same layer depth
input <- c(Sand[[x]] , Silt[[x]], Clay[[x]], OrgC[[x]], BulkDens[[x]], CEC[[x]], pHH2O[[x]], CrsFrg[[x]],ExchBases[[x]], EC[[x]], ExchAcids[[x]])
x = 1
input <- c(Sand[[x]] , Silt[[x]], Clay[[x]], OrgC[[x]], BulkDens[[x]], CEC[[x]], pHH2O[[x]], CrsFrg[[x]],ExchBases[[x]], EC[[x]], ExchAcids[[x]])
input
names(input) <- c("Sand", "Silt", "Clay",   "OrgC",  "BulkDens", "CEC" , "pHH2O", "CrsFrg", "ExchBases" , "EC", "ExchAcids")
input
VMCPTF <- list(NULL)
View(VMCPTF.function)
?exp
VMCPTF[[x]]
print(paste0("Done running ", "VMCPTF[[",x, "]]"))
VMCPTF <- list(NULL) # create an empty list to store the results per interval depth
for (x in 1:6) { # loop to apply the function to 6 depths
# compile in a raster stack all the input data needed for the same layer depth
input <- c(Sand[[x]] , Silt[[x]], Clay[[x]], OrgC[[x]], BulkDens[[x]], CEC[[x]], pHH2O[[x]], CrsFrg[[x]],ExchBases[[x]], EC[[x]], ExchAcids[[x]])
# adapt the names of the raster compiled layers
names(input) <- c("Sand", "Silt", "Clay",   "OrgC",  "BulkDens", "CEC" , "pHH2O", "CrsFrg", "ExchBases" , "EC", "ExchAcids")
# apply the function VMCPTF.function
VMCPTF[[x]] <- VMCPTF.function( Sand=input[["Sand"]],   Silt=input[["Silt"]], Clay=input[["Clay"]],
OrgC=input[["OrgC"]], BulkDens=input[["BulkDens"]], CEC=input[["CEC" ]],
pHH2O=input[["pHH2O"]],  h1=-10, h2=-20, h3=-31.6, h4=-1585,
PTF.coef, fix.values=TRUE, print.coef=TRUE)
print(paste0("Done running ", "VMCPTF[[",x, "]]"))
} # h0 is defined and then VMC_h0 is calculated in the function