library("sp")
library("ncdf")
library("raster")

# =======================
#-- get Gridded ERSST
# =======================

#- get Gridded ERSSTv3
# datadir = "../../data/ersst/data/ftproot/nc/ersst"
# file <- NULL
# i=1
# for (y in 1854:2010) {
# for (m in 1:12) {
#    file[i] = paste(datadir,y*100+m,"nc",sep=".")
#    i = i+1
# }}
# ERSST<- stack(file)

#- anomalize - 190101 - 200012
# start <- (1961-1854)*12+1
# end <- (1991-1854)*12

# ERSST_mask <- raster(ERSST,layer=1) # use for masking
# x <- getValues(ERSST_mask)
# x[x< -9] = NA # NAvalue = -9.99
# ERSST_mask = setValues(ERSST_mask,x)

# ERSST_bs = stack(unstack(ERSST)[start:end]) # get 30 years of 12 months
# bs = rep(1:12,30)   # mark each month
# ERSST_bs <- stackApply(ERSST_bs, bs, mean, na.rm=T) # get mean for each month in baseline
# ERSST <- ERSST - ERSST_bs # subtract mean from each month to make anomaly record
# rm (ERSST_bs)
# ersst <- mask(ERSST,ERSST_mask) # mask off the original NA 
# rm (ERSST)
# rm (ERSST_mask)

# ys <- as.vector(t(matrix(rep(seq(1,200),12),ncol=12)))
#- truncated to length of SST
# ys <- ys[1:dim(ersst)[3]]

#- annual means
# ersst <- stackApply(ersst,ys,mean)
# eryears <- ys[1:dim(ersst)[3]] 

#- write file for reuse
# writeRaster(ersst, filename='ersst-annual-anom-1961-1990.nc', format="CDF", overwrite=TRUE)

#- read prepared file (see above)
ersst <- brick("../../data/ersst/data/ersst-annual-anom-1961-1990.nc")


# =======================
#-- get Gridded HadSST2
# =======================

#- read the nc file
# hadfile = "../../data/hadcru/data/HadSST2_1850on.nc"

#- select sea surface temperatures
# HADSST<-brick(hadfile,varname='sst')

#-------------------------------------------
#-- Convert to Gridded Annual Anomalies
#-------------------------------------------

#ys <- as.vector(t(matrix(rep(seq(1,200),12),ncol=12)))
#- truncated to length of SST
# hadyears <- ys[1:dim(HADSST)[3]]
#- annual means
# hadsst <- stackApply(HADSST,hadyears,mean,na.rm=T,overwrite=TRUE)
# hadyears <- ys[1:dim(hadsst)[3]] 
# rm(HADSST)

# writeRaster(ersst, filename="hadsst-annual-anom-1961-2000.nc", format="CDF", overwrite=TRUE)
hadsst<-brick("../../data/hadcru/data/hadsst-annual-anom-1961-2000.nc")

# =======================

name = "NW Pacific"
region = "nwpac";
coo <- "23N-65N,120E-180E"
y1= 65; y2= 23; x1= 120; x2= 180; l1= 120; l2= 180 

# =======================

cropRaster <- function (r,r1,r2,c1,c2) {
        #r <- raster(r)
        return(crop(r, extent(r,rowFromY(r,r1), rowFromY(r,r2), colFromX(r,c1), colFromX(r,c2))))
}

region_hadsst <- stack(cropRaster(raster(hadsst,layer=1),y1,y2,x1,x2))
for (i in 2:dim(hadsst)[3]) {
    region_hadsst <- addLayer(region_hadsst,cropRaster(raster(hadsst,layer=i),y1,y2,x1,x2))

}

region_ersst <- stack(cropRaster(raster(ersst,layer=1),y1,y2,l1,l2))
for (i in 2:dim(ersst)[3]) {
    region_ersst <- addLayer(region_ersst,cropRaster(raster(ersst,layer=i),y1,y2,l1,l2))

}

# =======================

getRasterLatMean <- function(r) {

    rows = dim(r)[1]
    cols = dim(r)[2]

    wcell <- area(r,na.rm=F) # cells

    m <- mask(wcell,r)
    r1 <- r*wcell
    
    return(sum(extract(r1,1:ncell(r1)),na.rm=T)/sum(extract(m,1:ncell(m)),na.rm=T))

}

ersst_trend <- NULL
for (j in 1:dim(region_ersst)[3]) {
    ersst_trend[j] <- getRasterLatMean(raster(region_ersst,layer=j))
}

hadsst_trend <- NULL
for (j in 1:dim(region_hadsst)[3]) {
    hadsst_trend[j] <- getRasterLatMean(raster(region_hadsst,layer=j))
}

# =======================

library("Cairo")
library("maps")
source("plot.ts.spectrum.R")
source("plot.morlet.R")

hadsst.colors=c("#2909d8", "#264dff", "#3fa0ff", "#72d9ff", 
                "#aaf8ff", "#e0ffff", "#ffffbf", "#ffe099", 
                "#ffad72", "#f86d5e", "#d82632", "#a50021" )

ticks=c(-10,-5,-3,-1,-0.5,-0.2,0,0.2,0.5,1,3,5,10)

#- HadSST REGION
img = paste("img/trend-hadsst-region-",region,"-fig1.png",sep="")
CairoPNG(img,width=600,height=400)
r = raster(region_hadsst, layer=(dim(region_hadsst)[3]-1))
v = getValues(r)
nc = ncol(r)
nr = nrow(r)
mtit = paste("HadSST2 2010",name,"\n",coo)
plot(r, main=mtit,col=hadsst.colors,axis.args=list(at=ticks,labels=ticks),horizontal=T,breaks=ticks)
grid(nx=nc,ny=nr,col="lightgray")
grid(nx=nc/2,ny=nr/2,col="darkgray",lty="solid")
map("world",add=T)
dev.off()

#- ERSST REGION
img = paste("img/trend-ersst-region-",region,"-fig1.png",sep="")
CairoPNG(img,width=600,height=400)
r = raster(region_ersst, layer=(dim(region_ersst)[3]))
nc = ncol(r)
nr = nrow(r)
mtit = paste("ERSSTv3B 2010",name,"\n",coo)
plot(r, main=mtit,col=hadsst.colors,axis.args=list(at=ticks,labels=ticks),horizontal=T,breaks=ticks)
grid(nx=nc,ny=nr,col="lightgray")
grid(nx=nc/2,ny=nr/2,col="darkgray",lty="solid")
map("world",add=T)
dev.off()

#- HadSST WORLD
img = paste("img/trend-hadsst-region-",region,"-fig2.png",sep="")
CairoPNG(img,width=600,height=400)
r = raster(hadsst, layer=(dim(region_hadsst)[3]-1))
e <- extent(r, rowFromY(r,y1+1),rowFromY(r,y2),colFromX(r,x1),colFromX(r,x2))
nc = ncol(r)
nr = nrow(r)
mtit = paste("HadSST2 2010",name)
plot(r, main=mtit,col=hadsst.colors,axis.args=list(at=ticks,labels=ticks),horizontal=T,breaks=ticks)
plot(e, add=TRUE, col='green', lwd=4)
grid(nx=nc,ny=nr,col="lightgray")
grid(nx=nc/2,ny=nr/2,col="darkgray",lty="solid")
map("world",add=T)
dev.off()

#- ERSST WORLD
# world map is -180 to 180, ersst is -1 to 359
img = paste("img/trend-ersst-region-",region,"-fig2.png",sep="")
CairoPNG(img,width=600,height=400)
r = raster(ersst, layer=(dim(region_ersst)[3]))
e <- extent(r, rowFromY(r,y1),rowFromY(r,y2),colFromX(r,l1),colFromX(r,l2))
nc = ncol(r)
nr = nrow(r)
mtit = paste("ERSSTv3b 2010",name)
plot(r, main=mtit,col=hadsst.colors,axis.args=list(at=ticks,labels=ticks),horizontal=T,breaks=ticks)
plot(e, add=TRUE, col='green', lwd=4)
grid(nx=nc/4,ny=nr/4,col="darkgray",lty="solid")
map("world",add=T)
dev.off()

#- HadSST + ERSST Regional Linear Trend
img = paste("img/trend-sst-region-",region,"-fig3.png",sep="")
CairoPNG(img,width=600,height=400)
y <- 1854:2010
mtit = paste("SST",name,"\n",coo)
plot(hadsst_trend[5:(length(hadsst_trend)-1)]~y,type="l",
    main=mtit,xlab="",ylab="" )
lines(ersst_trend~y,col="red")
legend(1880,0.8,
    c("HadSST2","ERSSTv3b"),
	pch=18,lty=1,col=c("black","red"))
grid()
dev.off()

# HadSST Fourier + Morlet
img = paste("img/trend-hadsst-region-",region,"-fig4.png",sep="")
CairoPNG(img,height=900,width=600)
hyears = 1850:2011
mtit = paste("HadSST2",name,"1850-2011")
plot.ts.spectrum(hadsst_trend,hyears,0.18,mtit)
dev.off()
sp <- spectrum(hadsst_trend[!is.na(hadsst_trend)], log="no", plot=F)
1/sp$freq[sp$spec >= (max(sp$spec) * 0.18)]
# [1] 162.000000  81.000000  54.000000  40.500000  27.000000  20.250000
# [7]   7.363636   6.000000   5.586207   3.521739   3.176471   3.115385
#[13]   2.892857   2.842105   2.655738   2.454545   2.417910   2.382353
#[19]   2.250000   2.219178   2.189189   2.076923

# ERSST Fourier + Morlet
img = paste("img/trend-ersst-region-",region,"-fig4.png",sep="")
CairoPNG(img,height=900,width=600)
ery = 1854:2010
mtit = paste("ERSSTv3b",name,"1854-2010")
plot.ts.spectrum(ersst_trend,ery,0.1,mtit)
dev.off()
sp <- spectrum(ersst_trend[!is.na(ersst_trend)], log="no", plot=F)
1/sp$freq[sp$spec >= (max(sp$spec) * 0.1)]
# [1] 160.000000  80.000000  53.333333  40.000000  26.666667  20.000000
# [7]  10.000000   9.411765   7.619048   6.666667   6.153846   5.714286
#[13]   5.517241   4.444444   4.000000   3.478261