#library("RNetCDF") library("sp") library("ncdf") library("raster") #-- Get linear ERSSTv3 #glsrc="ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.ocean.90S.90N.df_1901-2000mean.dat" #nhsrc="ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.ocean.00N.90N.df_1901-2000mean.dat" #shsrc="ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.ocean.90S.00N.df_1901-2000mean.dat" # ocean SST_nhsrc="../../data/ncdc/data/monthly.ocean.00N.90N.df_1901-2000mean.dat" SST_shsrc="../../data/ncdc/data/monthly.ocean.90S.00N.df_1901-2000mean.dat" SST_glsrc="../../data/ncdc/data/monthly.ocean.90S.90N.df_1901-2000mean.dat" SST_fgl<-read.table(SST_glsrc,fill=TRUE,na.strings="-999.0000") SST_fnh<-read.table(SST_nhsrc,fill=TRUE,na.strings="-999.0000") SST_fsh<-read.table(SST_shsrc,fill=TRUE,na.strings="-999.0000") # annual means SST_gl <- aggregate(SST_fgl[,3],list(SST_fgl[,1]), FUN=mean, na.rm=T) SST_nh <- aggregate(SST_fnh[,3],list(SST_fnh[,1]), FUN=mean, na.rm=T) SST_sh <- aggregate(SST_fsh[,3],list(SST_fsh[,1]), FUN=mean, na.rm=T) rm(SST_fgl) rm(SST_fnh) rm(SST_fsh) #-- get Gridded ERSSTv3 # read the nc file 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 }} SST<- stack(file) start <- (1901-1854)*12+1 end <- (2001-1854)*12 SST_mask <- raster(SST,layer=1) # use for masking x <- getValues(SST_mask) x[x< -9] = NA # NAvalue = -9.99 SST_mask = setValues(SST_mask,x) SST_bs = stack(unstack(SST)[start:end]) # get 30 years of 12 months bs = rep(1:12,30) # mark each month SST_bs <- stackApply(SST_bs, bs, mean, na.rm=T) # get mean for each month in baseline SST <- SST - SST_bs # subtract mean from each month to make anomaly record SST <- mask(SST,SST_mask) # mask off the original NA rm (SST_bs) rm (SST_mask) # print for posterity # SST rows = dim(SST)[1] cols = dim(SST)[2] months = dim(SST)[3] # Global mask rgl <- raster(SST,layer=1) rgl <- setValues(rgl, c(rep(1,cols*rows))) # NH mask rnh <- raster(SST,layer=1) rnh <- setValues(rnh, c(rep(1,cols*rows/2),rep(NA,cols*rows/2))) # SH mask rsh <- raster(SST,layer=1) rsh <- setValues(rsh,c(rep(NA,cols*rows/2),rep(1,cols*rows/2))) #----------------------------------- #-- Function: Mean across Latitudes #----------------------------------- # add mask getRasterLatMean <- function(r) { rows = dim(r)[1] cols = dim(r)[2] wcell <- area(r,na.rm=F) # cells wlat <- extract(wcell,seq(1,((rows-1)*cols+1),cols)) # lat bands m <- mask(wcell,r) r1 <- r*wcell r2 <- extract(r1,ncell(r1)) return(sum(extract(r1,1:ncell(r1)),na.rm=T)/sum(extract(m,1:ncell(m)),na.rm=T)) #latm <- NULL #latm <- apply(matrix(getValues(r,1,nrows=rows),nrow=rows),1,mean,na.rm=T) #return(sum(latm*wlat,na.rm=T)/sum(wlat[!is.na(latm)])) #return(sum(latm)/length(latm)) } ############################################### # Convert Monthly Anomaly to Monthly Mean # Convert Monthly Mean to Annual Global Mean ############################################### #-- convert from absolute to anomalies #-- use individual months #-- convert to 1961-1990 to match HadSSTv2 #----------------------------------- #-- Get Monthly Mean across Latitudes #----------------------------------- glmm <- NULL nhmm <- NULL shmm <- NULL for (j in 1:dim(SST)[3]) { glmm[j] <- getRasterLatMean(mask(raster(SST,layer=j),rgl)) nhmm[j] <- getRasterLatMean(mask(raster(SST,layer=j),rnh)) shmm[j] <- getRasterLatMean(mask(raster(SST,layer=j),rsh)) } #------------------------------------------- #-- Convert Monthly Mean to Annual #------------------------------------------- gl_annmean1 <- apply(t(matrix(glmm,nrow=12)),1,mean,na.rm=T) nh_annmean1 <- apply(t(matrix(nhmm,nrow=12)),1,mean,na.rm=T) sh_annmean1 <- apply(t(matrix(shmm,nrow=12)),1,mean,na.rm=T) ############################################### # Convert Monthly Anomaly to Annual Anomaly # Convert Annual Anomaly to Annual Global Mean ############################################### #------------------------------------------- #-- Convert to Annual Anomalies #------------------------------------------- ys <- as.vector(t(matrix(rep(seq(1,200),12),ncol=12))) # truncated to length of SST ys <- ys[1:dim(SST)[3]] # annual means s <- stackApply(SST,ys,mean) #rm(SST) rm(ys) #------------------------------------------- #-- Convert to Annual Means #------------------------------------------- # for each year in the brick gl_annmean2 <- NULL nh_annmean2 <- NULL sh_annmean2 <- NULL for (j in 1:dim(s)[3]) { # calculate mean of each latitude band r <- raster(s,layer=j) gl_annmean2[j] <- getRasterLatMean(mask(r,rgl)) nh_annmean2[j] <- getRasterLatMean(mask(r,rnh)) sh_annmean2[j] <- getRasterLatMean(mask(r,rsh)) } rm(r) #------------------------------------------- # test fit gl, fit to 1880 to 2011 #------------------------------------------- gl_1 <- gl_annmean1[27:157] gl_2 <- gl_annmean2[27:157] SST_gl <- SST_gl[1:131,2] gl_3 = (gl_1 + gl_2)/2 sum((gl_3-SST_gl)^2) sum((gl_2-SST_gl)^2) sum((gl_1-SST_gl)^2) #------------------------------------------- # test fit nh, exclude 2011 #------------------------------------------- nh_1 <- nh_annmean1[27:157] nh_2 <- nh_annmean2[27:157] SST_nh <- SST_nh[1:131,2] nh_3 = (nh_1 + nh_2)/2 sum((nh_3-SST_nh)^2) sum((nh_2-SST_nh)^2) sum((nh_1-SST_nh)^2) #------------------------------------------- # test fit sh, exclude 2011 #------------------------------------------- sh_1 <- sh_annmean1[27:157] sh_2 <- sh_annmean2[27:157] SST_sh <- SST_sh[1:131,2] sh_3 = (sh_1 + sh_2)/2 sum((sh_3-SST_sh)^2) sum((sh_2-SST_sh)^2) sum((sh_1-SST_sh)^2) ############################################### # Plot Both Methods agains the Published Trend ############################################### library("Cairo") source("plot.ts.spectrum.R") source("plot.morlet.R") years = 1880:2010 CairoPNG("img/trend-ersst-global-fig1.R",height=400,width=600) plot(SST_gl~ years, type="l", main="ERSSTv3b Global Gridded -v- Linear",xlab="",ylab="",ylim=c(-0.8,0.8)) lines(years,gl_2,col="blue") lines(years,gl_1,col="red") lines(years,gl_2-gl_1,col="gray") legend(1880,0.8,c("Published Linear ersst2", "Gridded Monthly Anomalies -> Monthly Means -> Annual Means", "Gridded Monthly Anomalies -> Annual Anomalies -> Annual Means", "Difference between methods"), pch=18,lty=1,col=c("black","red","blue","gray")) grid() dev.off() CairoPNG("img/trend-ersst-global-fig2.R",height=400,width=600) plot(SST_nh~ years, type="l", main="ERSSTv3b Northern Gridded -v- Linear",xlab="",ylab="",ylim=c(-0.8,0.8)) lines(years,nh_2,col="blue") lines(years,nh_1,col="red") lines(years,nh_2-nh_1,col="gray") legend(1880,0.8,c("Published Linear ersst2", "Gridded Monthly Anomalies -> Monthly Means -> Annual Means", "Gridded Monthly Anomalies -> Annual Anomalies -> Annual Means", "Difference between methods"), pch=18,lty=1,col=c("black","red","blue","gray")) grid() dev.off() CairoPNG("img/trend-ersst-global-fig3.R",height=400,width=600) plot(SST_sh~ years, type="l", main="ERSSTv3b Southern Gridded -v- Linear",xlab="",ylab="",ylim=c(-0.8,0.8)) lines(years,sh_2,col="blue") lines(years,sh_1,col="red") lines(years,sh_2-sh_1,col="gray") legend(1880,0.8,c("Published Linear ersst2", "Gridded Monthly Anomalies -> Monthly Means -> Annual Means", "Gridded Monthly Anomalies -> Annual Anomalies -> Annual Means", "Difference between methods"), pch=18,lty=1,col=c("black","red","blue","gray")) grid() dev.off() #- annual means years=1854:2010 CairoPNG("img/trend-ersst-global-fig4.R",height=900,width=600) plot.ts.spectrum(gl_annmean1[1:157],years,"ERSSTv3b Global Gridded 1854-2010") dev.off() CairoPNG("img/trend-ersst-global-fig5.R",height=900,width=600) plot.ts.spectrum(nh_annmean1[1:157],years,"ERSSTv3b Northern Hemisphere Gridded 1854-2010") dev.off() CairoPNG("img/trend-ersst-global-fig6.R",height=900,width=600) plot.ts.spectrum(sh_annmean1[1:157],years,"ERSSTv3b Southern Hemisphere Gridded 1854-2010") dev.off() #- monthly means months=1:length(glmm) CairoPNG("img/trend-ersst-global-fig7.R",height=900,width=600) plot.ts.spectrum(glmm,months,"ERSSTv3b Global Gridded 1854-2010") dev.off() CairoPNG("img/trend-ersst-global-fig8.R",height=900,width=600) plot.ts.spectrum(nhmm,months,"ERSSTv3b Northern Hemisphere Gridded 1854-2010") dev.off() CairoPNG("img/trend-ersst-global-fig9.R",height=900,width=600) plot.ts.spectrum(shmm,months,"ERSSTv3b Southern Hemisphere Gridded 1854-2010") dev.off() quit()