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 = "Tropical West Pacific" region = "twpac" coo <- "23S-S5N,120E-180E" y1= 23; y2= -23; x1= 120; x2= 180; l1= 120; l2= 180 hyears = 1850:2011 # 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)))) } # rounding errors for crop? 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 twpac #hadsst_trend <- hadsst_trend[22:length(hadsst_trend)-1] # twpac, 1870-2010 to eliminate NaNs #ersst_trend <- ersst_trend[17:length(ersst_trend)] # twpac, 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 2011",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 2011",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 twpac 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(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") mtit = paste("HadSST2",name,hyr[1],"-",hyr[length(hyr)]) #hyr = 1870:2010 hyr = 1850:2011 plot.ts.spectrum(hadsst_trend,hyr,0.15,mtit) dev.off() sp <- spectrum(hadsst_trend[!is.na(hadsst_trend)], log="no", plot=F) 1/sp$freq[sp$spec >= 0.15] #[1] 162.000000 54.000000 32.400000 14.727273 9.000000 8.100000 6.230769 # [8] 5.785714 # ERSST Fourier + Morlet img = paste("img/trend-ersst-region-",region,"-fig4.png",sep="") CairoPNG(img,height=900,width=600) #ery = 1870:2010 ery = 1854:2010 mtit = paste("ERSSTv3b",name,ery[1],"-",ery[length(ery)]) plot.ts.spectrum(ersst_trend,ery,0.08,mtit) dev.off() sp <- spectrum(ersst_trend[!is.na(ersst_trend)], log="no", plot=F) 1/sp$freq[sp$spec >= 0.08] #[1] 160.000000 80.000000 53.333333 40.000000 20.000000 5.714286