A quick look at temperature anomaly distributions

R code to look at changing temp distributions – follow up to Hansen/Sato/Ruedy

Story submitted by commenter Nullius in Verba

There has been a lot of commentary recently on the new Hansen op ed, with associated paper. Like many, I was unimpressed at the lack of context, and the selective presentation of statistics. To some extent it was just a restatement of what we all already knew – that according to the record the temperature has risen over the past 50 years – but with extreme value statistics picked out to give a more alarming impression.

But beyond that, it did remind me of an interesting question I’ve considered before but never previously followed up, which was to ask how the distribution had actually changed over time. We know the mean has gone up, but what about the spread? The upper and lower bounds?

So I plotted it out. And having done so, I thought it might be good to share the means to do so, in case anyone else felt motivated to take it further.

I’m not going to try to draw any grand conclusions, or comment further on Hansen. (There will be a few mild observations later.) This isn’t any sort of grand refutation. Other people will do that anyway. For the purposes of this discussion, the data is what it is. I don’t propose to take any of this too seriously.

I’ll also say that I make no guarantees that I’ve done this exactly right. I did it quickly, just for fun, and my code is certainly not as efficient or elegant as it could be. If anyone wants to offer improvements or corrections, feel free.

—

I picked the HadCRUT3 dataset to look at partly because it doesn’t use the extrapolation that GISSTEMP does, only showing temperatures in a gridcell if there are actual thermometers there. But it was also because I’d looked at it before and already knew how to read it!

For various reasons it’s still less than ideal. It still averages things up over a month and a 5×5 degree lat/long gridcell. That loses a lot of detail and narrows the variances. From the point of view of studying heatwaves, you can’t tell if it was 3 C warmer for a month or 12 C warmer for a week. So it clearly doesn’t answer the question, but we’ll have a look anyway.

Then I picked the interval from 1900-1950 to define a baseline distribution. I picked this particular interval as a compromise – because the quality of the earliest data is very questionable, there being very few thermometers in most parts of the world, and because the mainstream claims often refer only to the post-1950 period as being attributable to man.

Of course, you have the code, so if you don’t like that choice you can pick a different period.

The first plot shows just the distribution for each month. Time is along the x-axis, temperature anomaly up the y-axis, and darker shading is more probable.

(From ‘HadCRUT3 T-anom dist 5 small.png’)

Because the outliers fade into invisibility, here is a plot of log-probability that emphasizes them more clearly.

(From ‘HadCRUT3 T-anom log-dist 20 small.png’)

You can see particularly from the second one that the incidence of extreme outliers hasn’t changed much. (But see later.)

You can also see that the spread is a lot bigger than the change. The rise in temperatures is perceptible, but still smaller than the background variation. It does look like the upper bound has shifted upwards – about 0.5 C by eye.

To look more precisely at the change, what I did next was to divide the distribution for each month by the baseline distribution for 1900-1950.

This says whether the probability has gone up or down. I then took logarithms, to convert the ratio to a linear scale. If doubling is so many units up, then halving will be the same number of units down. And then I colored the distribution with red/yellow for and increase in probability and blue/cyan for a decrease in probability.

(From ‘HadCRUT3 T-anom dist log-change 10 small.png’)

You can see now why it was convenient for Hansen to have picked 1950-1980 to compare against!

Blue surrounded by red means a broader distribution, as in the 1850-1870 period (small gridcell sample sizes give larger variances). Red surrounded by blue means a narrower distribution, as in the 1950-1990 period. Blue over red means cooler climate, as in the 1900-1920 period. Red over blue means a warmer climate, as in the post-1998 period.

The period 1900-1950 (baseline) shows a simple shift. The period 1950-1990 appears to be a narrowing and shift. The period post 1990 shows effects of shift meeting effects of narrowing. The reduction in cold weather since 1930 unambiguous, but increase in warm weather is only clear post 1998; until then, there was a decrease in warmer weather due to narrowing of distribution. There is a step change in 1998, little change thereafter.

Generally, the distribution has got narrower over the 20th century, but jumped to be wider again around the end of the century. (This may, as Tamino suggested, be because different parts of the world warm at different rates.) The assumption that the variance is naturally constant and any change in it must be anthropogenic is no more valid than the same assumption about the mean.

So far, there’s nothing visible to justify Hansen’s claims. The step-change after 1998 is a little bigger and longer than the 1940s, but not enough to be making hyperbolic claims of orders-of-magnitude increases in probability. So what have I missed?

The answer, it turns out, is that Hansen’s main results are just for summer temperatures. The spread of temperatures over land is vastly larger in the winter than it is the summer.

Looking at the summer temperatures only, the background spread is much narrower and the shifted distribution now pokes it head out above the noise.

I’ve shown the global figures for July below. I really ought to do a split-mix of July in the northern hemisphere and January on the southern, but given the preponderance of land (and thermometers!) in the NH, just plotting July shows the effect nicely.

(From ‘HadCRUT3 T-anom dist log-change 10 July.png’)

Again, the general pattern of a warm 1940s and a narrowing of the distribution from 1950-1980 shows up, but now the post-1998 step-change is now +2 C above the background. It also ramps up earlier, with more very large excursions (beyond 5 C) showing up around 1970, and a shift in the core distribution around 1988. The transitions look quite sharp. I think this is what Hansen is talking about.

Having had a quick look at some maps of where the hot spots are (run the code), it would appear a lot of these ‘heatwaves’ are in Siberia or northern Canada. I can’t see the locals being too upset about that… That also fits with the observation that a lot of the GISSTEMP global warming is due to the extrapolation over the Arctic. That would benefit from some further investigation.

None of this gives us anything solid about heatwaves, or tells us anything about cause, of course. For all we know, the same thing might have happened in the MWP.

Run the code! It generates a lot more graphics, at full size.

# ##################################################

 # R Script

 # Comparison of HadCRUT3 T-anom distributions to 1900-1950 average

 # NiV 11/8/2012

 # ##################################################

library(ncdf)

 library(maps)

# Local file location - ***CHANGE THIS AS APROPRIATE***

 setwd("C:/Data/Climate")

# Download file if no local copy exists

 if(file.exists("hadcrut3.nc") == FALSE) {

 download("http://www.metoffice.gov.uk/hadobs/hadcrut3/data/HadCRUT3.nc","hadcrut3.nc")

 }

# HadCRUT3 monthly mean anomaly dataset. For each 5x5 lat/long gridcell

 # reports the average temperature anomaly of all the stations for a

 # given month. Runs from 1850 to date, indexed by month.

 # i.e. month 1 is Jan 1850, month 25 is Feb 1852, etc.

 # Four dimensions of array are longitude, latitude, unknown, and month.

hadcrut.nc = open.ncdf("hadcrut3.nc")

# --------------------

 # Functions to extract distributions from data

# Names of months

 month = c("January", "February", "March", "April", "May", "June",

 "July", "August", "September", "October", "November", "December")

month_to_date = function(m) {

 mn = (m %% 12) + 1

 yr = 1850 + ((m-1) %/% 12)

 return(paste(month[mn]," ",yr,sep=""))

 }

# Function to show 1 month's data

 plotmonthmap = function(m) {

 d = get.var.ncdf(hadcrut.nc,"temp",start=c(1,1,1,m),count=c(72,36,1,1))

 clrs = rev(rainbow(403))

 brks = c(-100,-200:200/20,100)

 image(c(0:71)*5-180,c(0:35)*5-90,d,col=clrs,breaks=brks,useRaster=TRUE,

 xlab="Longitude",ylab="Latitude",

 main=paste("Temperature anomalies",month_to_date(m)))

 map("world",add=TRUE)

 }

# Function to extract one month's data, as a vector of length 36*72=2592

 getmonth = function(m) {

 res=get.var.ncdf(hadcrut.nc,"temp",start=c(1,1,1,m),count=c(72,36,1,1))

 dim(res) = NULL # Flatten array into vector by deleting dimensions

 return(res)

 }

# Given a vector of month indexes, extract data for all those months as a single vector

 getmultimonths = function(mvec) {

 res=vapply(mvec,getmonth,FUN.VALUE=c(1:2592)*1.0);dim(res)=NULL

 return(res)

 }

# Function to determine the T-anom distribution for a vector of month indexes

 # Result is a vector of length 402 representing frequency of 0.1 C bins

 # ranging from -20 C to +20 C. Element 200 contains the zero-anomaly.

 # Data is smoothed slightly to mitigate poor sample size out in the tails.

 gettadist = function(mvec) {

 d = getmultimonths(mvec)

 res = table(cut(d,c(-Inf,seq(-20,20,0.1),Inf)))[]

 res = res/sum(res,na.rem=TRUE)

 res = filter(res,c(1,4,6,4,1)/16)

 return(res)

 }

# --------------------

# Draw Summer maps 1998-2012

 for(m in c(seq(1781,1949,12),seq(1782,1949,12),seq(1783,1949,12))) {

 png(paste("HadCRUT3 T-anom map ",month_to_date(m),".png",sep="") ,width=800,height=500)

 plotmonthmap(m)

 dev.off()

 }

# Calculate average distribution 1900-1950

 # Late enough to have decent data, early enough to be before global warming

 baseline = gettadist(c(600:1199))

# Plot the distribution to see that it looks sensible

 png("HadCRUT3 T-anom 1900-1950 dist.png",width=800,height=500)

 plot(c(100:300)*0.1-20.1,baseline[100:300],type="l",

 xlab="Temperature Anomaly C",ylab="Frequency /0.1 C",

 main="Temperature Anomaly Distribution 1900-1950")

 dev.off()

# --------------------

# A few functions for plotting

# Add a semi-transparent grid to a plot

 gr = function() {abline(h=c(-20:20),v=seq(1850,2010,10),col=rgb(0.6,0.6,1,0.5))}

# Plot some data

 plotdist = function(data,trange=20,isyears=FALSE,dogrid=FALSE,main) {

 if(isyears) { drange = c(1:dim(data)[1]) + 1850 } else { drange = c(1:dim(data)[1])/12 + 1850 }

 trange = c((200-trange*10):(200+trange*10))

 image(drange,trange*0.1-20.1,

 data[,trange],

 col=gray(rev(0:20)/20),useRaster=TRUE,

 xlab="Year",ylab="Temperature Anomaly C",main=main)

 if(dogrid) { gr() }

 }

# Colour scheme and breaks for change plot

 # Scale is logarithmic, so goes from 1/7.4 to 1/2.6 to 1 to 2.7 to 7.4

 redblue=c(rgb(0,1,1),rainbow(10,start=3/6,end=4/6),rainbow(10,start=0,end=1/6),rgb(1,1,0))

 e2breaks=c(-100,-10:10/5,100)

# This generates a scale for the change plots

 png("E2 Scale.png",width=200,height=616)

 image(y=exp(c(-10:10/4)),z=matrix(c(-10:10/4),nrow=1,ncol=21),

 col=redblue,breaks=e2breaks,log="y",xaxt="n",ylab="Probability Ratio")

 dev.off()

plotdistchange = function(data,baseline,trange=20,isyears=FALSE,dogrid=FALSE,main) {

 if(isyears) { drange = c(1:dim(data)[1]) + 1850 } else { drange = c(1:dim(data)[1])/12 + 1850 }

 trange = c((200-trange*10):(200+trange*10))

 bdata=data/matrix(rep(baseline,each=dim(data)[1]),nrow=dim(data)[1])

 image(drange,trange*0.1-20.1,

 log(bdata[,trange]),

 col=redblue,breaks=e2breaks,useRaster=TRUE,

 xlab="Year",ylab="Temperature Anomaly C",main=main)

 if(dogrid) { gr() }

 }

# --------------------

# Make an array of every month's distribution

 datalst = aperm(vapply(c(1:1949),gettadist,FUN.VALUE=c(1:402)*1.0))

# Plot it out for a quick look

 for(w in c(3,5,10,20)) {

 png(paste("HadCRUT3 T-anom dist ",w,".png",sep=""),width=2000,height=600)

 plotdist(datalst,w,dogrid=TRUE,main="Temperature Anomaly Distribution")

 dev.off()

 }

# and a small one

 png("HadCRUT3 T-anom dist 5 small.png",width=600,height=500)

 plotdist(datalst,5,dogrid=TRUE,main="Temperature Anomaly Distribution")

 dev.off()

# Log-probability is shown to emphasise outliers

 png("HadCRUT3 T-anom log-dist 20.png",width=2000,height=600)

 plotdist(log(datalst),20,dogrid=TRUE,main="Temperature Anomaly Distribution (log)")

 dev.off()

# and a small one

 png("HadCRUT3 T-anom log-dist 20 small.png",width=600,height=500)

 plotdist(log(datalst),20,dogrid=TRUE,main="Temperature Anomaly Distribution (log)")

 dev.off()

# Now plot the change

png("HadCRUT3 T-anom dist log-change 20.png",width=2000,height=600)

 plotdistchange(datalst,baseline,20,dogrid=TRUE,main="Change in Temperature Anomaly Distribution")

 dev.off()

# Plot the middle +/-10 C, red means more common than 1900-1950 average, blue less common

 png("HadCRUT3 T-anom dist log-change 10.png",width=2000,height=600)

 plotdistchange(datalst,baseline,10,dogrid=TRUE,main="Change in Temperature Anomaly Distribution")

 dev.off()

# and a small one

 png("HadCRUT3 T-anom dist log-change 10 small.png",width=600,height=500)

 plotdistchange(datalst,baseline,10,dogrid=TRUE,main="Change in Temperature Anomaly Distribution")

 dev.off()

# --------------------------------

 # Analysis by month

# Reserve space

 datam = rep(0,162*402*12);dim(datam)=c(12,162,402)

 basem = rep(0,402*12);dim(basem)=c(12,402)

# Generate an array of results for each month, and plot distributions and changes

 for(m in 1:12) {

 datam[m,,] = aperm(vapply(seq(m,m+161*12,12),gettadist,FUN.VALUE=c(1:402)*1.0))

 basem[m,] = gettadist(seq(m+600,m+1200,12))

png(paste("HadCRUT3 T-anom dist 5 ",month[m],".png",sep=""),width=800,height=600)

 plotdist(datam[m,,],5,isyears=TRUE,dogrid=TRUE,

 main=paste(month[m],"Temperature Anomaly Distribution"))

 dev.off()

png(paste("HadCRUT3 T-anom log-dist 20 ",month[m],".png",sep=""),width=800,height=600)

 plotdist(log(datam[m,,]),20,isyears=TRUE,dogrid=TRUE,

 main=paste(month[m],"Temperature Anomaly Distribution (log)"))

 dev.off()

for(w in c(10,15,20)) {

 png(paste("HadCRUT3 T-anom dist log-change ",w," ",month[m],".png",sep=""),width=800,height=600)

 plotdistchange(datam[m,,],basem[m,],w,isyears=TRUE,dogrid=TRUE,

 main=paste(month[m],"Change in Temperature Anomaly Distribution"))

 dev.off()

 }

 }

# Done!

# --------------------------------

# Interpretation

 # --------------

 # Blue surrounded by red means a broader distribution, as in the 1850-1870 period

 # (small gridcell samples give larger variances).

 # Red surrounded by blue means a narrower distribution, as in the 1950-1990 period.

 # Blue over red means cooler climate, as in the 1900-1920 period.

 # Red over blue means a warmer climate, as in the post-1998 period.

# Observations

 # ------------

 # Period 1900-1950 (baseline) shows a simple shift.

 # Period 1950-1990 appears to be a narrowing and shift.

 # Period post 1990 shows effects of shift meeting effects of narrowing.

 # Reduction in cold weather since 1930 unambiguous, but increase in warm weather only clear post 1998;

 # until then, there was a decrease in warmer weather due to narrowing of distribution.

 # Step change in 1998, little change thereafter.

 # Post 2000 change is a bit bigger and longer lasting but not all that dissimilar to 1940s.

# Monthly

 # -------

 # Picture looks very different when split out by month!

 # Spread in summer is far narrower than in winter.

 # Offset in summer in 21st C exceeds upper edge of distribution from 1900-1950.

 # Excess over 1940s looks about 0.5 C.

 # Looking at maps, a lot of it appears to be in Siberia and Northern Canada.

 # There is still a narrowing of the distribution 1950-1990.

 # Step change occurs a little earlier, around 1990, to 1940s levels

 # then jumps again in 1998 to higher level. This is within 20th C annual bounds

 # but outside narrower summer bounds.

# Caveats

 # -------

 # The data is averaged over 5x5 degree gridcells monthly. Spread of

 # daily data would be 2-5 times broader (?).

 # Spread at point locations would be broader still.

 # Note also 'great dying of the thermometers' around 1990/2003 - could affect variance.

 # Has been reported prior to Soviet collapse cold weather exaggerated, to get fuel subsidy (?).

 # Distribution extremes poorly sampled, results unreliable beyond +/-5 C.

 # (May be able to do better with more smoothing?)

 # Temp anomalies outside +/-10 C occur fairly often, 'heatwaves' probably in this region.

 # No area weighting - gridcells nearer poles represent smaller area.

 # (This is intended to look at point distributions, not global average, though.

 # Area weighting arguably not appropriate.)

 # No error estimates have been calculated.

 # All the usual stuff about UHI, adjustments, homogenisation, LIA/MWP, etc.

# ##################################################