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

library(ncdf)

ausgabe <- c("./")

# ------------

#NCRUCELL <- 67420
#NCRUCELL <- 90881 # FMS / land_lad, including Antarctica and everything
#NCRUCELL <- 92113 # FMS / land_lad, including Antarctica and everything

zz <- file("grid.bin","rb")
# skip LPJ grid file header 35 Bytes
#seek(zz,where=35,origin="start")
#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) # longitude 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

# latitudes run -90 to 90
# longitudes run 0 to 180 and -180 to 0

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

#for(i in 1:NCRUCELL) {
#      print(paste("i ", i, " lon ", x[i*2-1], " ilon ", ilon_cru[i], " lat ", x[i*2], " ilat ", ilat_cru[i]))
#}

# ------------

# GLWD level 3 contains for each grid cell a type code
# in range 1..12; -1 == missing value
# Resolution is 30sec == 43200 x 21600 cells
# so that 60x60 cells of GLWD are aggregated into one 0.5x0.5deg cell for LPJ

frac.lake <- array(data=0,dim=c(720,360))
frac.res <- array(data=0,dim=c(720,360))
frac.river <- array(data=0,dim=c(720,360))
frac.marsh <- array(data=0,dim=c(720,360))
frac.swamp <- array(data=0,dim=c(720,360))
frac.coastal <- array(data=0,dim=c(720,360))
frac.saline <- array(data=0,dim=c(720,360))
frac.bog <- array(data=0,dim=c(720,360))
frac.interm <- array(data=0,dim=c(720,360))
frac.wet75 <- array(data=0,dim=c(720,360))
frac.wet37 <- array(data=0,dim=c(720,360))
frac.wet12 <- array(data=0,dim=c(720,360))

#nc.glwd <- open.ncdf("/p/projects/waterforce/heinke/GreenBlueBook/arrows/GLWD/glwd.nc")
nc.glwd <- open.ncdf("glwd_3.nc")

# latititudes run -90 to 90
# longitudes run -180 to 180

for(ilat in 1:360)
{
  print(paste("Process latitude band ",ilat,sep=""))
  chunk.glwd <- get.var.ncdf(nc.glwd,start=c(1,1+(ilat-1)*60),count=c(43200,60))

  for(ilon in 1:720)
  {
    cell <- chunk.glwd[(1:60)+(ilon-1)*60,]

    # meaning of data values taken from legend in glwd_3.avl    
    frac.lake[ilon,ilat] <- length(which(cell==1))/3600
    frac.res[ilon,ilat] <- length(which(cell==2))/3600
    frac.river[ilon,ilat] <- length(which(cell==3))/3600
    frac.marsh[ilon,ilat] <- length(which(cell==4))/3600
    frac.swamp[ilon,ilat] <- length(which(cell==5))/3600
    frac.coastal[ilon,ilat] <- length(which(cell==6))/3600
    frac.saline[ilon,ilat] <- length(which(cell==7))/3600
    frac.bog[ilon,ilat] <- length(which(cell==8))/3600
    frac.interm[ilon,ilat] <- length(which(cell==9))/3600
    frac.wet75[ilon,ilat] <- length(which(cell==10))/3600
    frac.wet37[ilon,ilat] <- length(which(cell==11))/3600
    frac.wet12[ilon,ilat] <- length(which(cell==12))/3600

  }
}

close.ncdf(nc.glwd)

frac.lake <- frac.lake[,360:1]
frac.res <- frac.res[,360:1]
frac.river <- frac.river[,360:1]
frac.marsh <- frac.marsh[,360:1]
frac.swamp <- frac.swamp[,360:1]
frac.coastal <- frac.coastal[,360:1]
frac.saline <- frac.saline[,360:1]
frac.bog <- frac.bog[,360:1]
frac.interm <- frac.interm[,360:1]
frac.wet75 <- frac.wet75[,360:1]
frac.wet37 <- frac.wet37[,360:1]
frac.wet12 <- frac.wet12[,360:1]

print(paste("dim(frac.lake) ", dim(frac.lake)))

# save.image("process_glwd.RData")


lakes <- array(dim=NCRUCELL)
reservoirs <- array(dim=NCRUCELL)
wetlands <- array(dim=NCRUCELL)
for(i in 1:NCRUCELL)
{
  #print(paste("i ", i, " lon ", x[i*2-1], " ilon_cru[i] ", ilon_cru[i] , " lat ", x[i*2], " ilat_cru[i] ", ilat_cru[i]))
  lakes[i] <- frac.lake[ilon_cru[i],ilat_cru[i]]
  #print(paste("i ", i, " lakes[i] ", lakes[i] ))
  reservoirs[i] <- frac.res[ilon_cru[i],ilat_cru[i]]
  wetlands[i] <- frac.wet12[ilon_cru[i],ilat_cru[i]]*0.125+frac.wet37[ilon_cru[i],ilat_cru[i]]*0.375+frac.wet75[ilon_cru[i],ilat_cru[i]]*0.75+frac.bog[ilon_cru[i],ilat_cru[i]]+frac.swamp[ilon_cru[i],ilat_cru[i]]+frac.marsh[ilon_cru[i],ilat_cru[i]]
}

#print("lakes array")
#print(lakes)
#print("lakes[which(lakes != 0)]")
#print(lakes[which(lakes != 0)])

# lake.bin is expected to contain only single bytes by the current LPJ trunk, fraction * 100
#writeBin(as.vector(lakes),paste(ausgabe,"glwd_lakes.bin",sep=""),size=4)
#writeBin(as.vector(reservoirs),paste(ausgabe,"glwd_reservoirs.bin",sep=""),size=4)
#writeBin(as.vector(wetlands),paste(ausgabe,"glwd_wetlands.bin",sep=""),size=4)

writeBin(as.vector(as.integer(lakes*100)),paste(ausgabe,"glwd_lakes.bin",sep=""),size=1)
writeBin(as.vector(as.integer(reservoirs*100)),paste(ausgabe,"glwd_reservoirs.bin",sep=""),size=1)
writeBin(as.vector(as.integer(wetlands*100)),paste(ausgabe,"glwd_wetlands.bin",sep=""),size=1)



colors <- rainbow(150,start=0,end=5/6)

library(fields)
#library(maptools)
#world_borders<-readShapePoly("/home/heinke/data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shp")

#bitmap(paste(ausgabe,"glwd_lakes.png",sep=""), height=10, width=20, res=200)
bitmap(paste(ausgabe,"glwd_lakes.png",sep=""), height=5, width=10, res=200)
#image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.lake, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.lake, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-90,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_reservoirs.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.res, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_river.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.river, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_marsh.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.marsh, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_swamp.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.swamp, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_coastal.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.coastal, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_saline.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.saline, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_bog.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.bog, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_interm.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.interm, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_wet50-100.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.wet75, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_wet25-50.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.wet37, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_wet0-25.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.wet12, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_wet.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.wet12*0.125+frac.wet37*0.375+frac.wet75*0.75, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()

bitmap(paste(ausgabe,"glwd_wetlands.png",sep=""), height=10, width=20, res=200)
image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), frac.wet12*0.125+frac.wet37*0.375+frac.wet75*0.75+frac.marsh+frac.swamp+frac.bog, zlim=c(1/3600,1), xlim=c(-180,180), ylim=c(-65,90), col=colors, nlevel=151, xlab="", ylab="")
#plot(world_borders,xlim=c(-180,180),ylim=c(-65,90),add=TRUE,fg="transparent")
dev.off()