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)
hadyears <- ys[1:dim(hadsst)[3]] 
rm(HADSST)

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

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

name = "NE Pacific"
region = "nepac"
coo <- "23N-65N,100W-180W"
y1= 65; y2= 23; x1= -180; x2= -100; l1= 180; l2= 260 
hyears = 1870:2010 # NaN at 1868

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

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))
}

# for NEPAC
hadsst_trend <- hadsst_trend[22:length(hadsst_trend)-1] # NEPAC, 1870-2010 to eliminate NaNs
ersst_trend <- ersst_trend[17:length(ersst_trend)] # NEPAC, just to align

plot(ersst_trend,type="l")
lines(hadsst_trend,col="red")
# =======================

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))
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="lightgray")
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
y <- 1870:2010 # for nepac
mtit = paste("SST",name,"\n",coo)
plot(hadsst_trend~y,type="l",
    main=mtit,xlab="",ylab="" )
lines(ersst_trend~y,col="red")
legend(1960,-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)
mtit = paste("HadSST2",name,"1870-2010")
hyr = 1870:2010
plot.ts.spectrum(hadsst_trend,hyr,0.1,mtit)
dev.off()
sp <- spectrum(hadsst_trend[!is.na(hadsst_trend)], log="no", plot=F)
 1/sp$freq[sp$spec >= 0.1]
# [1] 144.000000  72.000000  48.000000  36.000000  24.000000  20.571429
# [7]  18.000000   6.000000   5.538462   4.800000

# ERSST Fourier + Morlet
img = paste("img/trend-ersst-region-",region,"-fig4.png",sep="")
CairoPNG(img,height=900,width=600)
ery = 1870:2010
mtit = paste("ERSSTv3b",name,"1870-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/sp$freq[sp$spec >= 0.1]
# [1] 144.000000  72.000000  48.000000  36.000000  20.571429   6.000000   5.538462