library(rgdal)
library(rworldxtra)
library(raster)
data(countriesHigh)
library(raadtools)
library(rgeos)
library(RColorBrewer)
library(graticule)
# files with model output
fs <- list.files("/home/shared/data/krillhabitat/grids", pattern = "grd", full.names = TRUE)

prj<- CRS("+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=-40")

# map data
w <- spTransform(countriesHigh, prj)
# getData is broken, temporarily proceeed with no south georgia and falklands coastlines
#sgs<- spTransform(getData("GADM", country = c("SGS"), level = 0), prj)
#flk<- spTransform(getData("GADM", country = c("FLK"), level = 0), prj)

w2 <- spTransform(coastmap("ant_coast10"), prj)

# Colours and periods for plotting
warm.col<- "#B2182B"
cool.col<- "#2166AC"
periods <- c("ave11yr","ipcc2012")
the.pal <- brewer.pal(9,"YlOrRd")[3:9]
the.pal2<- paste0(the.pal, 40)

ex <- extent(-65,-22, -75,-55)
gridextent<- extent(-60,-25, -80,-50)

gridextent<- extent(projectExtent(raster(gridextent, crs = "+proj=longlat"), prj))
dummy <- as(raster(gridextent, nrow = 1200, ncol =1200, crs = prj), "SpatialPoints")

ex2<- projectExtent(raster(ex, crs = "+proj=longlat"), prj)
values(ex2)<- NA

op <- par(mfcol=c(1,2), mar=c(0,0,0,0), oma=c(1,2,2,1))
for(i in 1:2){ # first level of loop is for time period: 1st current then future
    hi.r <- brick(grep(paste(periods[i], "hi", sep = "_"), fs, value = TRUE))[[3]]
    sw.r <- brick(grep(paste(periods[i], "fswth", sep = "_"), fs, value = TRUE))[[3]]
    ai.r <-  brick(grep(paste(periods[i], "aice", sep = "_"), fs, value = TRUE))[[3]]
    rdg2.r <- brick(grep(paste(periods[i], "rdg2", sep = "_"), fs, value = TRUE))[[3]]

    # set filters: hi.r is thickness, sw.r is light availability at the bottom of the column; ai.ris concentration
    h <- (hi.r >= 1 & hi.r <= 2 & ai.r > 0.15) * rdg2.r # layer for ridging
    h <- projectRaster(from = h, crs = prj)

    h[!h > 0] <- NA_real_
    brks2<- c(0, 0.025, 0.05, 0.1, 0.2, 0.4, 0.8,1)

    plot(ex2, axes=F, asp = 1)

    ih <- as.image.SpatialGridDataFrame(as(h, "SpatialGridDataFrame"))
    image(ih, add = TRUE, col = the.pal, breaks = brks2)

    plot(w2, border = "gray60", col = c("lightcyan", "grey90")[w2$cst10srf], add = TRUE, lwd=0.5)
    #plot(sgs, add=T, col="gray90", border="gray60", lwd=0.25)
    #plot(flk, add=T, col="gray90", border="gray60", lwd=0.25)

    llgridlines(dummy, side = "EN", ndiscr = 80, plotLabels = T,
                easts = c(-70,-60, -50, -40, -30, -20, -10), norths = c(-80, -75,-70, -65, -60, -55), col="gray60", lwd=0.5, cex=0)

    box(col="gray60")

    #line for maximum ice extent
    ie <- rasterToContour(ai.r, level=0.15)
    plot(spTransform(ie, CRS(projection(prj))), add = TRUE, lty="dashed", col="lightblue2", lwd=2)

    # panel labelling
    if (i==1 ) mtext("Current", side=3, line =1, xpd=NA, col=cool.col)
    if (i==2 ) mtext("Warm", side=3, line =1, xpd=NA, col=warm.col)
        if (i==1) mtext("September", side=2, las=0,outer=F, line=1, xpd=NA)

    # manual colourbar
    # increment y position of fill by 0.2e6; put text at bottom of each
    if (i==2){
      op2<- par(xpd=NA)
      for(j in 1:length(the.pal2)){

      inc<- 0.15e6
      adj<- (j-1)*inc
      xl<- 1.3e6
      xr<-xl+1e5
      yb<- 2.5e+06+adj
      yt<- yb+0.15e6
      tcx<- 0.4
      ofs=0.2

      rect(xleft=xl, ybottom=yb, xright=xr, ytop=yt, col= the.pal[j], border=NA)
      text(xr,yb, brks2[j], pos=4, cex=tcx, offset=ofs)
      if(j==(length(the.pal)))text(xr,yb+inc, brks2[j+1], pos=4, cex=tcx, offset=ofs)
      }
      par(op2)
    }
  }
par(op)