I did this work several weeks ago. It’s the same as Roman’s hammer except that there is no offset. Statistically this is inferior to the method that Roman presented, however this is the standard in climatology. Nick Stokes pointed out several times that I needed to do the anomaly correctly for this post, I chose the period 1950-1980. The code to do the weighting is slightly more sophisticated than previous posts and more accurate. Previously, 5×5 degree polar gridcells were weighted by area and averaged for the whole earth. The weights were normalized such that missing cells didn’t affect trend. In this version each year is weighted by grid and availability so that years with different patterns of missing data are appropriately weighted for trend.
This post is an apples to apples comparison to the methods presented in Roman’s hammer. Same data, same gridding, different (standard but less accurate) anomaly combination methods.
rgl=gridtem #rgl = global temperature
dim(rgl)=c(36*72,2532) #reset dimensions
rgl=t(rgl) #transpose for weighting and processing
rgl=ts(rgl,start=1800,deltat=1/12) #make into time series
mask=!is.na(colMeans(rgl,na.rm=TRUE)) #make mask for removing gridcells with no data
weight=sin(seq(2.5,180, 5)*(pi)/180) #cell weights by latitude
wt=rep(weight,72) #copy weights across whole longitude of globe
maskb=array(FALSE,dim=c(36,72)) #mask for hemisphere
maskb[1:18,]=TRUE #set sh first
dim(maskb)=(36*72) #dimension mask for combination
maskc=mask & maskb #combine masks
wt=wt/mean(wt[maskc]) #scale weights according to hemisphere's available data
fullmask=is.na(rgl) #set mask for missing data
fullwt=rep(wt,2532) #repeat global grid weight mask for all years -- R uses a ton of memory
dim(fullwt)=(c(2592,2532)) #set up full array for weights
fullwt=t(fullwt) #transpose array to match temp matrix - done for clarity
fullwt[fullmask]=NA #set full matrix missing values to NA
meanwt=rowMeans(fullwt,na.rm=TRUE) #take mean of each row for renormilization of weights for missing values per year
fullwtnorm=fullwt[1:2532,]/meanwt #renormalize
sh=ts(rowMeans(rgl[,maskc] * fullwtnorm[,maskc],na.rm=TRUE),start=1800,deltat=1/12) #calc southern hemisphere
maskc=mask & !maskb #flip mask to northern hemisphere
nh=ts(rowMeans(rgl[,maskc] * fullwtnorm[,maskc],na.rm=TRUE),start=1800,deltat=1/12) #calc northern hemisphere
un=ts.union(nh,sh) #unionize time series
glt=ts(rowMeans(un,na.rm=TRUE),start=1800,deltat=1/12) #take mean for global temperature
If you don’t regularly run R code, it won’t read easily but the comments give the steps. Global annual trend is below.
Read the rest of this entry »
