# R program to plot both SST and Ice coverage on the same figure

library(grid)
library(PP.map)

pp.map.internal.getPalette.icescale<-colorRampPalette(
       c(
         rgb( 255, 255, 255 ,   0, maxColorValue = 255 ),
         rgb( 255, 255, 255 , 255, maxColorValue = 255 )
        )
) 

f=pp.ppa('/ibackup/cr1/hadobs/OBS/marine/HadISST/anoms/HadISST1.1_sst_1870on_5dg_anm6190.pp',
         max=1)
colours=pp.map.internal.getPalette.diverging(10)
getGp<-function(x) {
    minV = min(f[[1]]@data,na.rm=T)
    maxV = max(f[[1]]@data,na.rm=T)
    fraction=(x-minV)/(maxV-minV)
    colour=colours[as.integer(10*fraction)+1]
    return(gpar(fill=colour,col=colour))
}
pp.map.internal.wm <- map('world',interior=FALSE,plot=FALSE)
is.na(pp.map.internal.wm$x[8836])=T  # Remove Antarctic bug

pushViewport(viewport(width=0.9,height=0.5,x=0.05,y=0.25,
                      just=c("left","bottom"),name="vp_map"))
pushViewport(plotViewport(margins=c(2,2,2,2)))
pushViewport(dataViewport(c(-180,180),c(-90,90)))
polygons = pp.get.poly(f[[1]])
#gp2= getGp(2.0)
for(i in seq(1,f[[1]]@lbrow)) {
   for(j in seq(1,f[[1]]@lbnpt)) {
      if(!is.na(f[[1]]@data[i*f[[1]]@lbnpt+j])) {
#         print(sprintf("%d %d\n",i,j))
#         print(polygons[1,,j,i])
#         print(polygons[2,,j,i])
         grid.polygon(x=polygons[1,,j,i],y=polygons[2,,j,i],
                      default.units="native",
                      gp=getGp(f[[1]]@data[i*f[[1]]@lbnpt+j]))
      }
   }
}
ice=pp.ppa('/ibackup/cr1/hadobs/OBS/marine/HadISST/ice/HadISST1_ice_1870on.pp',
         max=1)
colours=pp.map.internal.getPalette.icescale(10)
getGp<-function(x) {
    minV = min(ice[[1]]@data,na.rm=T)
    maxV = max(ice[[1]]@data,na.rm=T)
    fraction=(x-minV)/(maxV-minV)
    colour=colours[as.integer(10*fraction)+1]
    return(gpar(fill=colour,col=colour))
}
polygons = pp.get.poly(ice[[1]])
for(i in seq(1,ice[[1]]@lbrow)) {
   for(j in seq(1,ice[[1]]@lbnpt)) {
      if(!is.na(ice[[1]]@data[i*ice[[1]]@lbnpt+j])) {
#         print(sprintf("%d %d\n",i,j))
#         print(polygons[1,,j,i])
#         print(polygons[2,,j,i])
         grid.polygon(x=polygons[1,,j,i],y=polygons[2,,j,i],
                      default.units="native",
                      gp=getGp(ice[[1]]@data[i*ice[[1]]@lbnpt+j]))
      }
   }
}
llines(pp.map.internal.wm$x,pp.map.internal.wm$y,col="black")
grid.xaxis(main=T)
grid.xaxis(main=F)
grid.yaxis(main=T)
grid.yaxis(main=F)

upViewport(0)

makeGPars <- function(n) {

}