So I’ve been trying to do a little more QC with the GHCN data. The problem with the raw data is that a single station ID has multiple curves and very little explanation as to what causes the differences. All of the following curves were taken from the same instrument. Yes they have different data, yes they have odd steps compared to each other and unfortunately careful examination reveals that the data is a mess.



If you take the anomaly of this data and then average the trend has substantially less slope than if you first average then take the anomaly. Since these are the same data, it makes sense to take the anomaly after averaging. If you have different datasets, it makes sense to take the anomaly before averaging.

I spent about an hour this morning wrestling with Phil Jone’s version of GHCN which is somewhat different. I was able to determine that he used in the case above an average first then anomaly combination, which in this case made sense.

To that end, I’ve redone my own gridded global temperature. It’s still not properly area weighted but each step gets closer. I focused primarily on the proper identification of sub-timeseries inside GHCN temperature stations which are actually from the same instrument.

Here is my improved GetStation algorithm for GHCN data compilation. You pass the right station id and it performs some QC to determine the correct method to combine the data.

#combine staton data by insuring whether each series is in fact the same

getstation4=function(staid=91720)

{

alladj=NA

allraw=NA

#raw data

mask= (ghmean[,2]==staid)

data=ghmean[mask,]

noser=levels(factor((data[,3])))

for(i in noser)

{

mask2=data[,3]==i

startyear=min(data[mask2,4])

endyear=max(data[mask2,4])

subd=data[mask2,5:16]

index=(data[mask2,4]-startyear)+1

dat=array(NA,dim=c(12,(endyear-startyear)+1))

for(j in 1:length(index))

{

dat[,index[j]]=as.numeric(subd[j,])

}

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))>=2010)

{

print (“error series”)

rawd=NA

} if(!is.ts(allraw))

{

allraw=rawd

}

else

{

allraw=ts.union(allraw,rawd)

}

} nc=ncol(allraw)

matchv=array(NA,dim=c(nc,nc))

mv=rep(NA,nc)

if(nc>1)

{

for (i in 1:nc)

{

for (j in 1:nc)

{

dt=allraw[,i]-allraw[,j]

mask=!is.na(dt)

if (sum(mask,na.rm=TRUE)>2)

{

if (mean(abs(dt),na.rm=TRUE)<.2)

{#data is the same

matchv[i,j]=j

}

}

}

}

index=1

for (i in 1:nc)

{

mask=!is.na(matchv[i,])

if( sum(!is.na(mv[mask]))>0 )

{#found match to already matched number

mv[mask]=min(mv[mask],na.rm=TRUE)

}else{

mv[mask]=index

index=index+1

}

}

}else{

mv=1

} srs=levels(factor(mv))

outar=array(NA,dim=c(nrow(allraw),length(srs)))

index=1 for (i in srs)

{

mask=i==mv

if(!is.null(ncol(allraw[,mask])))

{

outar[,index]=rowMeans(allraw[,mask],na.rm=TRUE)

}else{

outar[,index]=allraw[,mask]

} index=index+1

}

outar=ts(outar,start=time(allraw)[1],deltat=1/12) outara=calc.anom(outar)

if(!is.null(ncol(outara)))

{

final=ts(rowMeans(outara,na.rm=TRUE),start=time(outara)[1],deltat=1/12)

}else{

final=outara

}

final

}

I hate what wordpress does to the formatting. The first section separates the data into timeseries as shown in Figure 1. This section was improved to recognize skipped years in the data.

After the data is in each series, they are all compared against each other to determine how similar they are. Since not many versions of the EXACT SAME data are an exact match, I used this line.

if (mean(abs(dt),na.rm=TRUE)<.2)

This line takes the mean of the absolute value of the difference. If the difference is on average greater than 0.2 it decides that it is not the same data. Consider what that means for our ability to detect 0.5C/century trends. I know of at least one instance where my eyes told me it was the same data yet it failed this test. I just wasn’t willing to loosen it further. Really it shouldn’t make much difference but wow it’s noisy stuff considering that it’s from the same instrument at the same time. Why is there any question at all whether the thermometer measured 25.5 or 25.1?

Anyway, after the algorithm decides that the data is the same, all of the ‘same’ data is averaged into a single time series before conversion to anomaly as you would sensibly do if you knew it was from the same instrument.

The last section of the algorithm is a bit complicated in it’s function. In the 8 series above, you get a list like series 4 is the same as series 1,3, 8 and in series 8 you get 1,4 so you need to combine evrything so that we know series 3 is also the same as 1,4 and 8. Anyway, when it’s all done if there are multiple “substantially different” series, they are then anomalized and combined into a single GHCN time series.

The data is then gridded and the north and south hemisphere are averaged separately.

A comparison to CRU is a little bit problematic as this doesn’t include any ocean data — 70% of the earth. It also doesn’t demonstrate the huge 1998 spike, don’t know why. Overall the pre-1970 data is substantially higher than the CRU version and consequently the trend has reduced century scale warming to about 0.45 C/century, although without ocean data this version shows more warming since 1980. It doesn’t make sense to look at century trend to me, but there has been less warming than advertised according to this curve.

At this point, there needs to be a lot more work on data QC and another look at urban vs rural. There are too many steps, oddities and station deletions in recent years to make any conclusions at this point but it seemed worthwhile to present the early results.



