Thermal Hammer

Ok, today is the day that tAV, the repeatedly alleged EVIL denialist blog, lucky unwitting recipient of the biggest scientific scandal of the last hundred years, hosting the proprietor who has been snipped from every advocate website, who will now present a global gridded temperature from the RAW-ish GHCN data, having a higher trend than the believers will publish.  However, it is to my knowledge a more correct representation of the actual data.  To perform this feat of magic, we’ll use Roman’s hammer, so popular it has been immortalized in a t-shirt.

Roman’s hammer, is a simple method for the combination of temperature time series. Remember different thermometers can experience two primary effects. Offset of temperature due to altitude or proximity to water etc.  and different levels of seasonal variance, proximity to water, dryness of air, etc.. The “hammer” method takes care of both seasonal variation and offsetting values to provide the best match between series. A second and independent improvement from Roman’s recent work means that we no longer need to calculate anomaly in order to solve for global temperature trend.

The steps are simple.

– load data
– sort temperature series into their own 5×5 gridcells 72lon by 36 lat
– create individual series from the GHCN inventory file by simple averaging of station ID’s with multiple series representing the same instrument.
– gather all series for each gridcell and hammer them together.
– average NH an SH individually and combine — not sure if this makes sense but it copies HadCrut.

I’ve changed the renown “getstation” algorithm for this post, it has a new value called “version”.  Since multiple instruments which sometimes have the same WMO number are mixed in with multiple copies of the same instrument (How sloppy is that?!) they use an additional digit. Considering the crew, team or whatever, wants literally trillions of dollars to limit CO2, they can certainly spend the small time required to document and record “actual” temperature stations by their own ID.   In this post, the algorithm sorts out which stations are nearby vs same instrument.  Multiple series are placed in columns and averaged using a row mean and placed in a single output timeseries.

I’ll put the whole code at a different link but here is the getsinglestation function:

getsinglestation=function(staid=60360, version=0)
{
	staraw=NA
	#raw data
	smask= (ghmean[,2]==staid)
	mask= ghmean[smask,3]==version
	data=ghmean[smask,][mask,]
	noser= levels(factor(data[,4]))

	for(j in noser)
	{
		mask2=data[,4]==j
		startyear=min(data[mask2,5])
		endyear=max(data[mask2,5])
		subd=data[mask2,6:17]
		index=(data[mask2,5]-startyear)+1
		dat=array(NA,dim=c(12,(endyear-startyear)+1))
		for(k in 1:length(index))
		{
			dat[,index[k]]=as.numeric(subd[k,])
		}
		dim(dat)=c(length(dat),1)

		dat[dat==-9999]=NA
		dat=dat/10
		rawd=ts(dat,start=startyear,deltat=1/12)
		if(max(time(rawd))>=2011)
		{
			print ("error series")
			rawd=NA
		}

		if(!is.ts(staraw))
		{
			staraw=rawd
		}else{
			staraw=ts.union(staraw,rawd)
		}
	}
	if(!is.null(ncol(staraw)))
	{
		allraw=ts(rowMeans(staraw,na.rm=TRUE),start=time(staraw)[1],freq=12)
	}else{
		allraw=staraw
	}

	allraw
}

The main loop of the software is worth discussing a bit. In the code below tinv is the inventory values for each temp station  It contains lat, lon, and a variety of station metadata. Before the loop,  an array is created to which gridded cell each station is assigned to.

##assign stations to gridcells
gridval=array(NA,dim=c(stationnum,2))
gridval[,1]=as.integer((tinv[,5]+90)/5)+1  #lat
gridval[,2]=as.integer((tinv[,6]+180)/5)+1  #lon

Over the last couple of months, the algorithm keeps getting simpler, which of course makes the engineer in me happier with each revision.   It starts where “ii” is gridded longitude in 5 degree increments and “jj” is latitude also each 5 degrees.

I’ve added a number of comments to the code for this post.

gridtem=array(NA,dim=c(36,72,211*12))

for(ii in 1:72)
{
	for(jj in 1:36)
	{
		maska=gridval[,2]==ii  #mask all stations where longitude doesn't match
		maskb=gridval[,1]==jj  #mask all stations where latitude doesn't match
		maskc= maska & maskb   #combine lon and lat masks to get all stations where both match
		sta=NA                 ##initialize to non TS
		for(i in (1:stationnum)[maskc])
		{
			rawsta = getsinglestation(tinv[i,2],tinv[i,3])

			if(!is.ts(sta)) #add stations to multicolumn timeseries
			{
				sta=rawsta            #first station
			}else{
				ts.union(sta,rawsta)  #other columns
			}
		}

		if(is.ts(sta))           #if at least one time series exists
		{

			sta=window(sta,start=1800)	                #trim any pre-1800 data
			index=as.integer(((time(sta)-1800)*12+.002)+1)  #calculate time index for insertion into array
			if(!is.null(ncol(sta)))                         #is station more than one column
			{
				gridtem[jj,ii,index]=temp.combine(sta)$temps    #more than one column, use Roman's algorithm
			}else{
				gridtem[jj,ii,index]=sta                        #a single column, just assign the data
			}
			print(jj)                                       #progress debug
			tit=paste("ii=",ii,"jj=",jj)
			plot(ts(gridtem[jj,ii,],start=1800,freq=12),main=tit)   #debug plot
		}
	}
	print ("COLUMN")        #progress debug
	print(ii) 		#progress debug
}

That’s it. Pretty simple, just gather it into gridcells and if there is more than one temp station per gridcell, use the Roman Hammer to smack them into place. Ok Jeff, we’ve waited for two weeks, enough scribbling of nonsense, what are the results.

Remember Roman’s temperature combinations allow different offsets for linear trend by month, yet a single trend for the whole dataset.  The offsets are required to align anomalies with each other and represent a substantial improvement over typical global instrumental temperature series.

The algorithm above is called with the following lines.

rgl=gridtem                                                #move gridded temperature to different variable
dim(rgl)=c(36*72,2532)                                     #restructure grid into 36*72 individual series
rgl=t(rgl)                                                 #transform rows and columns so each year is a row
rgl=ts(rgl,start=1800,deltat=1/12)                         #make time series
mask=!is.na(colMeans(rgl,na.rm=TRUE))                      #create mask for empty columns (empty gridcells)
wt=rep(weight,72)                                          #apply area weighting to gridcells
maskb=array(FALSE,dim=c(36,72))
maskb[1:18,]=TRUE                                          #create mask by hemisphere
dim(maskb)=(36*72)
maskc=mask & maskb                                         #combine hemisphere mask with empty gridcell mask
nh=temp.combine(rgl[,maskc],wt[maskc])$temps               #use Roman's method to combine all gridcell series in the hemisphere
maskc=mask & !maskb                                        #invert hemisphere mask and redo for sh
sh=temp.combine(rgl[,maskc],wt[maskc])$temps

glt=(nh+sh)/2                                              #take global average

It’s pretty sweet how simple things are getting. There is of course still the possibility (an in my case probability) of error but as the code has been gone over again and again, the differences created by any mistakes are becoming small.

Northern hemisphere trend.

Figure 1 - Northern hemisphere temperatures

I bet you’ve never seen a semi-global temperature plot that looks like that!  It happens to be in Degrees J,  which have equal increments to degrees C but the true value has some offset. Think of it as Degrees C, except that nobody knows where zero really is.

The Southern hemisphere is in Figure 2.

Figure 2 - Southern hemisphere temperatures

Pretty unique looking. It’s a simple matter to combine the two series, yup I used the hammer again.

glt=temp.combine(un)$temps

Figure 3, Global temperatures

Ok, so that told us nothing, except that we’ve now been able to calculate global temperature in degrees J.  The varition in actual value seems fairly small.

Trend is what we care about tho.

Using Roman’s true anomaly methods, I get the following plot.

Figure 4, Global temperature trend with true slope

This method really is different from the less accurate yet standard anomaly trend of climate science.  Globally, 0.079 C/Decade since 1900.  Of course this is just GHCN data, so the trend could be created by UHI and such.

Below is a corrected anomaly trend.

Figure 4 Global temperature trend by anomaly

Below that is the same plot as Fig 4 with HadCrut overlaid.The anomaly for this post is calculated using my adaptation of Roman’s trend for anomaly calculation.

Figure 5 - Global temperatures plus HadCrut

Now look at that match.  CRU has slightly too low a historic value WRT  GHCN 1935 and before, but in the alleged anthropogenic global warming era, GHCN has a much higher trend.  I call it Id’s but it’s really primarily Roman’s work which has been shaped into this post.  There are better methods for area combination than gridded and we should play with those in time, but there is a lot of information in the above graph. As a last minute add on, Bob Tisdale has a post which is a fantastic match to this result.

First the obvious, a skeptic, denialist, anti-science blog published a greater trend than Phil Climategate Jones.  What IS up with that?

From climategate, we learned that Phil flinched away from review of at least one technical math oriented paper.  My own guess is that he just hadn’t considered an offset method for aligning anomalies just hoping the steps would come out in the wash.   In reality, when there is an upslope in the true signal, the steps created from starting and stopping of anomaly temperature series, always reduce the trend.

We discussed that in this post Anomaly Aversion.

The hemispheres tell an interesting story.

Figure 6 - Northern hemisphere temp anomaly

Figure 6 - Southern hemisphere temperature anomaly

Check out the difference between the two hemispheres of GHCN.

Figure 7 - Both hemispheres. Filtered 11 year Gaussian

Note the delta in the most recent data.

There are a lot of details in this next plot which others may not see.  Many of you are more experienced and knowledgeable with respect to temp trend than I, but running the algorithm makes a difference in understanding the quality of the data.  Visually, as each dataset is plotted by anomaly, I see a general upslope in most data.  There are plenty of blue series which go against the trend though, which brings the question of data quality into focus.

Below is a global plot since 1978.  The deep red and blue points are gridcells with such unreasonable temperature trends that they cannot be accepted.  There are far more deep red than deep blue.  So far though, this algorithm is still a bit of a black box beast.  I’ve not explored the reasons for individual extremes in gridcells and we are looking for small temp trends.

Figure 8 - Temperature trends since 1978

Figure 9 - Temperature since 1900

This post has been worked on for literally weeks, I’ve looked at so many plots my head spins.  There is a lot of unwritten detail here.  Consider the amazing evenness in Russian temp stations.  Enough so that it will be difficult to ever question Siberia as was done prior to climategate.

There are high trends from GHCN,  so high in fact that anyone who questions Phil Climategate Jones temp trends  will need to show some evidence. Certainly Phil is an ass, but it no longer seems to me that he has ‘directly’ exaggerated temp trends one bit.   Also, the elimination of temp stations has certainly reduced the quality of data, however, despite E.M. Smiths work, we don’t know which way the bias runs.   We need the actual data for that.

Several skeptics will dislike this post.  They are wrong, in my humble opinion.   While winning the public “policy” battle outright,  places pressure for a simple unified message, the data is the data and the math is the math.  We”re stuck with it, and this result.  In my opinion, it  is a better method.  Remember though, nothing in this post discusses the quality of the raw data. I’ve got a lot of information on data quality, for the coming days.  In the meantime, consider what would cause such a huge difference in trend between the northern and southern hemispheres.

Anyway, the global temp trend, from this data since 1978 is:

Figue 10, global temperatures as calculated from GHCN