convert texture
# load required libraries
library(soiltexture)
library(caTools)
library(ncdf)
in_directory<-"../Soil_BNU/CN/nc/"
silt=read.ENVI(paste(in_directory,"ENVI/SI",sep=""))
#clay=read.ENVI(paste(in_directory,"ENVI/CL",sep=""))
#sand=read.ENVI(paste(in_directory,"ENVI/SA",sep=""))
rows<-nrow(silt)
cols<-ncol(silt)
layers<-8
my_data <- data.frame("CLAY" = as.vector(clay[,,]),"SILT" = as.vector(silt[,,]),"SAND" = as.vector(sand[,,]))
summary(my_data)
my_data2<-my_data
rm(clay,silt,sand)
mask<-which(my_data[,1]>=0 & my_data[,2]>=0 & my_data[,3]>=0 & (my_data[,1]+my_data[,2]+my_data[,3])>50)
# Nomalized data to make the sum of Clay, Silt and Sand equate to 100
my_data2[mask,] <- TT.normalise.sum(tri.data = my_data[mask,])
summary(my_data2)
code<-my_data2[,1]
code[-mask]<-NA
# get texture code
code[mask]<-TT.points.in.classes(tri.data = my_data2[mask,],class.sys = "USDA.TT", PiC.type = "t") #
# double class only the first class were used
code[grep(",",code)]<-unlist(strsplit(code[grep(",",code)],','))[c(seq(1, length(unlist(strsplit(code[grep(",",code)],','))), by=2))]
summary(as.factor(code))
texture<-array(code, c(rows,cols,layers))
texture_code<-array(0, c(rows,cols,layers))
texture_code[texture=="Cl"]<-1
texture_code[texture=="SiCl"]<-2
texture_code[texture=="SaCl"]<-3
texture_code[texture=="ClLo"]<-4
texture_code[texture=="SiClLo"]<-5
texture_code[texture=="SaClLo"]<-6
texture_code[texture=="Lo"]<-7
texture_code[texture=="SiLo"]<-8
texture_code[texture=="SaLo"]<-9
texture_code[texture=="Si"]<-10
texture_code[texture=="LoSa"]<-11
texture_code[texture=="Sa"]<-12
texture_code_t<-aperm(texture_code, c(2,1,3))
texture_code_t<-texture_code_t[4320:1,,]
write.ENVI(texture_code_t,"texture_code_t")
#write.ENVI(texture_code,"texture_code")
# this is used for calculate top soil depth of the global soil database
# the thickness for each layer is (45,46,75,123,204,336,554,913)
# The start for each layer is (0,45,91,166,289,493,829,1383)
library(caTools)
library(lattice)
in_directory<-"../Soil_BNU/CN/Soil Hydraulic Parameters/nc/ENVI/"
# read required original data from "Soil Hydraulic Parameters"
As=read.ENVI(paste(in_directory,"THSCH",sep = ""))
Afld=read.ENVI(paste(in_directory,"TH33",sep = ""))
Awlt=read.ENVI(paste(in_directory,"TH1500",sep = ""))
Cs=read.ENVI(paste(in_directory,"PSI_S",sep = ""))
Ks=read.ENVI(paste(in_directory,"K_SCH",sep = ""))
# read required original data from "Soil Database"
in_directory<-"F:/Ning/WaSSI/Soil_BNU/CN/nc/out/"
texture=read.ENVI(paste(in_directory,"texture_code",sep = ""))
Depth_max=read.ENVI(paste(in_directory,"PDEP1",sep = ""))
texture<-texture[,,1:7] # only first 7 layers were used
# try to use depth data
Depth_max<-Depth_max*10 #Transfer depth form cm to mm
Depth_max[Depth_max<0]<-0
Depth_max[is.na(Depth_max)]<-0
summary(as.vector(Depth_max),na.rm=T)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 0.0 357.7 828.2 1702.0
# get dimision
rows<-nrow(texture)
cols<-ncol(texture)
layers<-7
# soil depth for each layer
layer_thickness<-c(45,46,75,123,204,336,554)
layer_start<-c(0,45,91,166,289,493,829,1383)
# this is the SINX= 5.08(1000/CN-10), Where Zup_depth= SINX/(As-Afld)
SINX<-array(0,c(rows,cols,layers))
SINX[texture<=3 & texture>0]<-134.5 # A HSG
SINX[texture<=6 & texture>=4]<-55.2 # B HSG
SINX[texture==7]<-30.2 # C HSG
SINX[texture<=12 & texture>=8]<-21.6 # D HSG
# Zup_depth
zup_depth<- SINX/(As-Afld)
# for each layer if the calculte soil depth > layer thickness, then set layer thickness as the default thickness
for (i in 1:3) {
zup_depth[,,i][zup_depth[,,i]>layer_thickness[i]]<-layer_thickness[i]
}
# since the fourth layer, if the texture is rock the set the layer depth to 0 for each layer
for (i in 4:7) {
zup_depth[,,i][texture[,,i]>8]<-0
zup_depth[,,i][zup_depth[,,i]>layer_thickness[i]]<-layer_thickness[i]
}
# base on Max soil depth to calculate the real up and sub depth of soil
zup<-matrix(0,nrow=rows,ncol=cols)
# Depth_max < Thickness to layer 1
class1<- which(Depth_max <= layer_start[2])
zup[class1]<- Depth_max[class1]/2
# Depth_max < Thickness to layer 2
class2<- which(Depth_max > layer_start[2] & Depth_max <= layer_start[3])
zup[class2]<-zup_depth[,,1][class2]
# Depth_max < Thickness layer 3
class3<- which(Depth_max > layer_start[3] & Depth_max <= layer_start[4])
zup[class3]<-(zup_depth[,,1]+zup_depth[,,2])[class3]
# Depth_max < Thickness layer 4
class4<- which(Depth_max > layer_start[4] & Depth_max <= layer_start[5])
zup[class4]<-(zup_depth[,,1]+zup_depth[,,2]+zup_depth[,,3])[class4]
# Depth_max < Thickness layer 5
class5<- (Depth_max > layer_start[5] & Depth_max <= layer_start[6])
zup[class5]<-(zup_depth[,,1]+zup_depth[,,2]+zup_depth[,,3]+zup_depth[,,4])[class5]
# Depth_max < Thickness layer 6
class6<- which(Depth_max > layer_start[6] & Depth_max <= layer_start[7])
zup[class6]<-(zup_depth[,,1]+zup_depth[,,2]+zup_depth[,,3]+zup_depth[,,4]+zup_depth[,,5])[class6]
# Depth_max < Thickness layer 7
class7<-which(Depth_max > layer_start[7] & Depth_max <= layer_start[8])
zup[class7]<-(zup_depth[,,1]+zup_depth[,,2]+zup_depth[,,3]+zup_depth[,,4]+zup_depth[,,5]+zup_depth[,,6])[class7]
# Depth_max > Thickness layer 7
class8<- which(Depth_max > layer_start[8])
zup[class8]<-(zup_depth[,,1]+zup_depth[,,2]+zup_depth[,,3]+zup_depth[,,4]+zup_depth[,,5]+zup_depth[,,6]+zup_depth[,,7])[class8]
# get low profile depth
zlow<-Depth_max-zup
#summary
summary(as.vector(zup))
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 0.0 0.0 207.6 459.3 1322.0 28089
summary(as.vector(zlow))
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 0.0 0.0 149.7 254.7 1453.0 28089
levelplot(zup)
levelplot(zlow)
levelplot(zup+zlow)
levelplot(Depth_max)
#write.ENVI(zup+zlow, "depth")
write.ENVI(zup, paste(in_directory,"zup",sep = ""))
write.ENVI(zlow, paste(in_directory,"zlow",sep = ""))