rm(list=ls(all=T))
gc()

library(fields)
library(soiltexture)

load("process_hwsd2.RData")

#soilclasses in the same order as in soil.par
soilclasses <- c("Cl","SiCl","SaCl","ClLo","SiClLo","SaClLo","Lo","SiLo","SaLo","Si","LoSa","Sa","ROCK")

#remove NAs
data.t_gravel[which(is.na(data.t_gravel))] <- 0
data.t_sand[which(is.na(data.t_sand))] <- 0
data.t_silt[which(is.na(data.t_silt))] <- 0
data.t_clay[which(is.na(data.t_clay))] <- 0

#NCRUCELL <- 67420
#NCRUCELL <- 90881 # FMS / land_lad, including Antarctica and everything
NCRUCELL <- 92113 # for CM2.1p1 coastlines

zz <- file("grid.bin","rb")   #CRUGRID
## skip LPJ grid file header 35 Bytes
#seek(zz,where=35,origin="start")
#x <- array(0,dim=2*NCRUCELL)
##x <- readBin(zz, integer(), n=2*NCRUCELL, size=2, endian="swap") / 100
#x <- readBin(zz, integer(), n=2*NCRUCELL, size=2) / 100

# read header
# headername version
# LPJGRID\3
hheadername <- readChar(zz, 7, TRUE)
hversion <- readBin(zz, integer(), n=1)
print(paste("Header: headername ", hheadername, " version ", hversion))
horder <- readBin(zz, integer(), n=1)     # order of data items CELLYEAR,YEARCELL or CELLINDEX
hfirstyear <- readBin(zz, integer(), n=1) # first year for data
hnyear <- readBin(zz, integer(), n=1)     # number of years
print(paste("Header: order ", horder, " firstyear ", hfirstyear, " nyear ", hnyear))
hfirstcell <- readBin(zz, integer(), n=1) # index of first data item
hncell <- readBin(zz, integer(), n=1)     # number of data item per year
hnbands <- readBin(zz, integer(), n=1)    # number of data elements per cell
hcellsize_lon <- readBin(zz, double(), n=1, size=4) # latitude cellsize in deg
hscalar <- readBin(zz, double(), n=1, size=4)    # conversion factor
print(paste("Header: firstcell ", hfirstcell, " ncell ", hncell, " nbands ", hnbands, " cellsize_lon ", hcellsize_lon, " scalar ", hscalar))

if (hversion == 3) {
   hcellsize_lat <- readBin(zz, double(), n=1, size=4) # latitude cellsize in deg
   hdatatype <- readBin(zz, integer(), n=1)
}

NCRUCELL <- hncell
x <- array(0,dim=hnbands*NCRUCELL)
x <- readBin(zz, integer(), n=hnbands*NCRUCELL, size=2) * hscalar
#print(paste(x))

cellarea <- (111e3*hcellsize_lon)*(111e3*hcellsize_lon)*cos(x[c(1:NCRUCELL)*2]/180*pi)
ilon <- as.integer((x[c(1:NCRUCELL)*2-1]+180)/hcellsize_lon + 1.01)
ilat <- as.integer((x[c(1:NCRUCELL)*2]+90)/hcellsize_lon + 1.01)
close(zz)

#for(i in 1:NCRUCELL) {
#      print(paste("i ", i, " ilon ", ilon[i], " ilat ", ilat[i]))
#}

#
vector.sea <- array(dim=NCRUCELL)
vector.soil <- array(dim=NCRUCELL)
vector.sand <- array(dim=NCRUCELL)
vector.silt <- array(dim=NCRUCELL)
vector.clay <- array(dim=NCRUCELL)
vector.ice <- array(dim=NCRUCELL)
vector.rock <- array(dim=NCRUCELL)
for(i in 1:NCRUCELL)
{
  vector.sea[i] <- frac.sea[ilon[i],ilat[i]]
  #merge "Dunsess and shifting sand" into soil assuming "sand" from Cosby, 1984
  vector.soil[i] <- frac.soil[ilon[i],ilat[i]] + frac.DS[ilon[i],ilat[i]]
  vector.sand[i] <- (data.t_sand[ilon[i],ilat[i]]*frac.soil[ilon[i],ilat[i]] + 92*frac.DS[ilon[i],ilat[i]])/(frac.soil[ilon[i],ilat[i]] + frac.DS[ilon[i],ilat[i]])
  vector.silt[i] <- (data.t_silt[ilon[i],ilat[i]]*frac.soil[ilon[i],ilat[i]] + 5*frac.DS[ilon[i],ilat[i]])/(frac.soil[ilon[i],ilat[i]] + frac.DS[ilon[i],ilat[i]])
  vector.clay[i] <- (data.t_clay[ilon[i],ilat[i]]*frac.soil[ilon[i],ilat[i]] + 3*frac.DS[ilon[i],ilat[i]])/(frac.soil[ilon[i],ilat[i]] + frac.DS[ilon[i],ilat[i]])
  vector.rock[i] <- frac.RK[ilon[i],ilat[i]]
  vector.ice[i] <- frac.GG[ilon[i],ilat[i]]	
  if (x[i*2] < -55.0){
    #print(i)
    vector.ice[i] <- 1
  }
}

#fill missing cell with global average particle size distribution
vector.sand[which(is.na(vector.sand))] <- sum(vector.sand*vector.soil,na.rm=T)/sum(vector.soil)
vector.silt[which(is.na(vector.silt))] <- sum(vector.silt*vector.soil,na.rm=T)/sum(vector.soil)
vector.clay[which(is.na(vector.clay))] <- sum(vector.clay*vector.soil,na.rm=T)/sum(vector.soil)


soil.data <- array(dim=c(NCRUCELL,3),dimnames=list(c(1:NCRUCELL),c("SAND","SILT","CLAY")))
soil.data[,1] <- vector.sand
soil.data[,2] <- vector.silt
soil.data[,3] <- vector.clay

soil.class <- TT.points.in.classes(tri.data=soil.data,class.sys="USDA.TT",PiC.type="l")

#soil.class[which(rowSums(soil.class)>1),]

#handle borderline cases
soil.class[which(soil.class[,"Lo"] & soil.class[,"SiLo"]),"SiLo"] <- FALSE
soil.class[which(soil.class[,"SaClLo"] & soil.class[,"SaLo"]),"SaLo"] <- FALSE
soil.class[which(soil.class[,"Cl"] & soil.class[,"ClLo"]),"Cl"] <- FALSE
soil.class[which(soil.class[,"SaLo"] & soil.class[,"LoSa"]),"SaLo"] <- FALSE

#image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), data.t_oc, zlim=c(0,max(data.t_oc,na.rm=T)), xlim=c(-180,180), ylim=c(-90,90), col=rainbow(n=50,start=0,end=5/6), xlab="",ylab="")


#zz <- file("soil.bin","rb")
#soil.old <- readBin(zz,integer(),n=NCRUCELL,size=1)
#close(zz)

soil.new <- array(data=0,dim=NCRUCELL)
for(i in 1:(length(soilclasses))-1) soil.new[which(soil.class[,soilclasses[i]])] <- i
soil.new[which(vector.rock>(vector.soil))] <- 13
#print(vector.ice)
soil.new[which(vector.ice>(vector.soil))] <- 14
#print(soil.new)

writeBin(as.vector(as.integer(soil.new)),"soil_new.bin",size=1)


#soil.new59199 <- soil.new
#soil.new59199[which(soil.old==0)] <- 0
#writeBin(as.vector(as.integer(soil.new59199)),"soil_new59199.bin",size=1)

#soil.new_skiprock <- soil.new
#soil.new_skiprock[which(soil.new_skiprock==13)] <- 0
#writeBin(as.vector(as.integer(soil.new_skiprock)),"soil_new_skiprock.bin",size=1)

## for visualisation:

zz <- file("soil_new.bin","rb")
soil.dummy <- readBin(zz,integer(),n=NCRUCELL,size=1)
close(zz)

soil.map <- array(dim=c(720,360))
for(i in 1:NCRUCELL) soil.map[ilon[i],ilat[i]] <- soil.dummy[i]

bitmap("soil_new.png", type="png16m", height=5, width=10, res=400, pointsize=10)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), soil.map, zlim=c(0,20), xlim=c(-180,180), ylim=c(-90,90), col=rainbow(n=50,start=0,end=5/6), xlab="",ylab="")
dev.off()