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

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
cellarea <- (111e3*0.5)*(111e3*0.5)*cos(x[c(1:NCRUCELL)*2]/180*pi)
ilon <- as.integer((x[c(1:NCRUCELL)*2-1]+180)/0.5 + 1.01)
ilat <- as.integer((x[c(1:NCRUCELL)*2]+90)/0.5 + 1.01)
close(zz)

#
#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)
vector.lake <- array(dim=NCRUCELL)
vector.csgl <- array(dim=NCRUCELL)

print(frac.WR)
print(frac.CSGL[which(frac.CSGL > 0)])

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
  #}
  vector.lake[i] <- frac.WR[ilon[i],ilat[i]] * 100
  vector.csgl[i] <- frac.CSGL[ilon[i],ilat[i]] * 100
}

print(vector.lake)
print(vector.csgl[which(vector.csgl > 0)])

##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

bitmap("frac.WR.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), frac.WR, zlim=c(0,max(frac.WR,na.rm=T)), xlim=c(-180,180), ylim=c(-90,90), col=rainbow(n=50,start=0,end=6/6), xlab="",ylab="")
dev.off()

bitmap("frac.CSGL.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), frac.CSGL, zlim=c(0,max(frac.CSGL,na.rm=T)), xlim=c(-180,180), ylim=c(-90,90), col=rainbow(n=50,start=0,end=6/6), xlab="",ylab="")
dev.off()

#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(vector.lake)),"lake_frac.bin",size=1)
writeBin(as.vector(as.integer(vector.csgl)),"great_lakes_frac.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("lake_frac.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("lake_frac.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,255), xlim=c(-180,180), ylim=c(-90,90), col=rainbow(n=50,start=0,end=5/6), xlab="",ylab="")
dev.off()