Early note on May temperatures

There is only a few hours left of May here, and the NCAR reanalysis has three days of data to go. May has had fewer extremes than earlier months, but warmed toward the end. It is currently showing the same average anomaly as February, and I think that is about where it will end up. That is actually not quite as warm as last May, or as March 2015, but it continues a warm start to the year. The GISS anomaly for February was 0.8°C.

GWPF submission

Last month, I posted (here and here) about a GWPF call for submissions on temperature adjustments. I included a hasty draft, which people were kind enough to comment on. I hope I've incorporated most of those suggestions in this update. I have edited for sobriety and brevity, and shortened the blogpost inclusions. I have added a list of references. I have left the in-line links, because I doubt that many readers will have a paper-only copy; there is a list of urls at the end.

I'm planning to submit it within the next few days, unless readers think further improvement is required.

 Update - I have sent it with the correction noted in comments.

Daily Arctic Temperatures

People are watching the Arctic Sea Ice, which is melting rather fast this year. Often quoted is a plot of temperatures above 80°N from DMI. I'll show it below.

Update New post here. maintained location here.

That plot is based on reanalysis data, and covers just one region. Since I am now regularly integrating current NCEP/NCAR reanalysis data (here), I thought I could do the same for that region, and the comparison would be a check. However, in some ways that high Arctic is not where melting is currently concentrated, so I thought I could also give a better spread of areas. Naturally, that means a gadget, below the fold.

The DMI plot shows an average from 1958-2002. That isn't without cost - it requires melding a number of different reanalyses to cover the range. It could be argued too that it misses recent Arctic warming, and so may show modern data as unusually warm, when by modern standards it isn't.

Anyway, I decided to stick with NCEP/NCAR data from 1994 to present, over which time it seems fairly reliable. And of course I use NCEP for current. There are some discrepancies. DMI aims to show the true surface temp 2m above the ground. I show what is basically the bottom grid level.

Here are the comparative plots for 2015. I've tried to stick to the DMI style here:

DMI	NCEP/NCAR vai Moyhu

I have however used Celsius - I can't see the point of Kelvin here, especially with freezing being prominent. The plots are broadly similar. The DMI average has a longer period above zero, but to a smaller extent. This may reflect the slight difference in levels. DMI shows a warmer looking winter - reflecting I think the more ancient reference average.

So below the fold is the gadget. It lets you choose arbitrary rectangles on a 5×15° grid in the region above 60°N. I'll place it just below, with further description below that.



The "rectangle" shows with colored sides. Below is a box of controls. The colored triangles move the sides in the indicated directions by one step. Outward pointing triangles expand the area. The black triangles translate the rectangle by longitude. The initial rectangle has overlapping sides and doesn't show all the colors.

When you have selected the area, press "Plot new" to show the plot of daily average temperatures for 2015 to date. Day of year is marked on the x axis.

I'll add the plot to the maintained latest data page.

How bad is naive temperature averaging?

In my last post, I described as "naive averaging" the idea of calculating an average temperature anomaly G by simply subtracting from each of a number of local records the (varying) lifetime average, and averaging the resulting differences. There, and in an earlier post I gave simple examples of why it didn't work. And in that last post, I showed how the naive average could be made right by iteration.

The underlying principle is that in making an anomaly you should subtract your best estimate of what the value would be. That leaves the question of how good does "best" have to be; it has to be good enough to resolve the thing you are trying to deduce. If that is the change of global temperature, your estimate has to be accurate to the effect of that change.

If you add a global G to a station mean, then the mean of the result isn't right unless mean G is zero. There is freedom to set average G over an interval, but only one - not over all station intervals. So as in the standard method, you can set G to have mean zero over a period like 1961-90, and use station means over that period as the offsets. Providing there are observations there, which is the rub. But much work has been done on methods for this, which itself is evidence that the faults of naive averaging are well known.

Anyway, here I want to check just how much difference the real variation of intervals in temperature datasets makes, and whether the iterations (and TempLS) do correct it properly. I take my usual collection of GHCN V3 and ERSST data, at 9875 locations total, and the associated area-based weights, which are zero for all missing values. But then I change the data to a uniformly rising value - in fact, equal to the time in century units since the start in 1899. That implies a uniform trend of 1 C/century. Do the various methods recover that?

Update - I made a small error in trends below. I used unweighted regressions, which allowed the final months of 2015 to be entered as zeroes. I have fixed this, and put the old tables in faint gray beside. The plots were unaffected. I have not corrected the small error I noted in the comments, but now it is even smaller - a converged trend of 0.9998 instead of 1. 

 So I'll start with land only (GHCN) data - 7280 stations - because most of the naive averaging is done with land stations. Here is the iterative sequence:

Iter Trend RMS 1 0.547165 4.079662 2 0.781293 1.23042 3 0.89165 0.529091 4 0.945758 0.251682 5 0.972713 0.12399 6 0.986226 0.061898 10 0.999021 0.003962 20 0.999899 4e-06

Iter Trend RMS 1 0.535192 4.07986 2 0.777397 1.230401 3 0.889767 0.529098 4 0.944707 0.251689 5 0.972055 0.123994 6 0.985764 0.0619 10 0.998745 0.003962 20 0.999636 4e-06

As in the last post, the first iteration is the naive calculation, where the offsets are just the means over the record length. RMS is a normalised sqrt sum squares of residuals. As you see, the trend for that first step was about half the final value. And it does converge to 1C/Cen, as it should. This is the value TempLS would return.

I should mention that I normalised the global result at each step relative to the year 2014. Normally this can be left to the end, and would be done for a range of years. Because of the uniformity of the data, setting to a one year anomaly base is adequate. For annual data, base setting is just the addition of a constant. For monthly data, there are twelve different constants, so the anomaly base setting has the real effect of adjusting months relative to each other. It makes a small difference within each year, but if not set till the end, the convergence of the trend to 1 is not so evident, though the final result is the same.

So normalised to 2014, here is the plot of the global temperatures at each iteration:



The first naive iteration deviates quite substantially, with reduced trend. The deviation is due solely to the pattern of missing values.

Now I'll do the same calculation including sea surface temperature - the normal TempLS range. The effects are subdued, because SST grid values don't generally have start and end years like met stations, even though they may have missing values.

Iter Trend RMS 1 0.828193 5.189208 2 0.968292 0.625476 3 0.993962 0.100924 4 0.998822 0.018172 5 0.999756 0.003425 6 0.999936 0.000658 10 0.999979 1e-06 20 0.999979 0

Iter Trend RMS 1 0.814107 5.189285 2 0.966456 0.625463 3 0.993393 0.100923 4 0.99844 0.018172 5 0.999406 0.003425 6 0.999592 0.000658 10 0.999637 1e-06 20 0.999637 0

The corresponding plot of iterative G curves is:



The deviation of the first step is reduced, but still considerable. Convergence is faster.

How to average temperature over space and time

This post is about the basics of TempLS. But I want to present it in terms of first principles of averaging, with reference to other methods that people try that just don't work. I'll mention the least squares approach described here.

I once described a very simple version of TempLS. In a way, what I am talking about here is even simpler. It incorporates an iterative approach which not only gives context to the more naive approaches, which are in fact the first steps, but is a highly efficient solution process.

Preliminaries - averaging methods

The first thing done in TempLS is to take a collection of monthly station readings, not gridded, and arrange them in a big array by station, month and year. I denote this x_smy. Missing values are marked with NA entries. Time slots are consecutive. "Station" can mean a grid cell for SST, assigned to the central location.

I'm going to average this over time and space. It will be a weighted average, with a weight array w_smy exactly like x. The first thing about w is that it has a zero entry wherever data is missing. Then the result is not affected by the value of x, which I usually set to 0 there for easier linear algebra in R.

The weighting will generally be equal over time. The idea is that averaging for a continuum (time or space) generally means integration. Time steps are equally spaced.

Spatial integration is more complicated. The basic idea is to divide the space into small pieces, and estimate each piece by local data. In effect, major indices do this by gridding. The estimate for each cell is the average of datapoints within, for a particular time. The grid is usually lat/lon, so the trig formula for area is needed. The big problem is cells with no data. Major indices use something equivalent to this.

Recently I've been using a finite element style method - based on triangular meshing. The surface is completely covered with triangles with a datapoint at each corner (picture here). The integral is that of the linear interpolation within each triangle. This is actually higher order accuracy - the 2D version of the trapezoidal method. The nett result is that the weighting of each datapoint is proportional to the area of the triangles of which it is a vertex (picture here).

Averaging principles

I've described the mechanics of averaging (integrating) over a continuum. Now it's time to think about when it will work. We are actually estimating by sampling, and the samples can be sparse.

The key issue is missing values. There is one type of familiar missing value - where a station does not report in some month. As an extension, you can think of the times before it started or after it closed as missing values.

But in fact all points on Earth except the sample points can be considered missing. I've described how the averaging deals with them. It interpolates and then integrates that. So the issue with averaging is whether that interpolation is good.

The answer for absolute temperatures is,no it isn't. They vary locally within any length scale that you can possibly sample. And anyway, we get little choice about sampling.

I've written a lot about where averaging goes wrong, as it often does. The key thing to remember is that if you don't prescribe infilling yourself, something else will happen by default. Usually it is equivalent to assigning to missing points the average value of the remainder. This may be very wrong.

A neat way of avoiding most difficulties is the use of anomalies. Here you subtract your best estimate of the data, based on prior knowledge. The difference, anomaly, then has zero expected value. You can average anomalies freely because the default, the average of the remainder, will also be zero. But it should be your best estimate - if not, then there will be a systematic deviation which can bias the average when missing data is defaulted. I'll come back to that.

What TempLS does

I'll take up the story where we have a data array x and weights w. As described here, TempLS represents the data by a linear model
x_smy=L_sm+G_my+ε_smy

Here L is a local offset for station month (like a climate normal), G is a global temperature function and ε an error term. Suffices show what they depend on.

TempLS minimises the weighted sum of squares of residuals:
Σw_smyd_smy² where d_smy=x_smy-L_sm-G_my

It is shown here that leads to two equations in L and G which TempLS solves simultaneously. But here I want to give them an intuitive explanation. The local offsets don't vary with year. So you can get an estimate by averaging the model over years. That gives an equation for each L.

And the G parameters don't involve the station s. So you can estimate those by averaging over stations. That would be a weighted average with w.

Of course it isn't as simple as that. Averaging over time to get L involves the unknown G. But it does mean, at least, that we have one equation for each unknown. And it does allow iteration, where we successively work out L and G, at each step using the latest information available. This is a well-known iterative method - Gauss-Seidel.

That's the approach I'll be describing here. It works well, but also is instructive. TempLS currently uses either block Gaussian elimination with a final direct solve of a ny×ny matrix (ny=number years), or conjugate gradients, which is an accelerated version of the iterative method.

Naive methods

A naive method that I have inveighed against is to simply average the absolute temperatures. That is so naive that it is, of course, never done with proper weighting. But weighting isn't the problem.

But one that is common, and sometimes not too bad, is where a normal is first calculated for each station/month over the period of data available, and used to create anomalies that are averaged. I gave a detailed example of how that goes wrong here. Hu McCulloch gave one here.

The basic problem is this. In setting T=L+G as a model, there is an ambiguity. L is fixed over years, G over space. But you could add any number to the L's as long as you subtract it from G. If you want L to be the mean over an interval, you'll want G to have zero mean. But If stations have a variety of intervals, G can't be zero-mean over all of them.

The normal remedy is to adopt a fixed interval, say 1961-90, on which G is constrained to have zero mean. Then the L's will indeed be the mean for each station for that interval, if they have one. The problem of what to do if not has been much studied. There's a discussion here.

Iterating

I'll introduce a notation Σw_s, Σw_y, meaning that a w-weighted sum is formed over the suffix variable (station, year) holding the others fixed. I'd call it an average, but the weights might not add to 1. So Σw_y yields an array with values for each station and month. Our equations are then
Σw_y(d_smy) = 0
or Σw_y(L_sm) = Σw_y(T_smy - G_my)......Eq(1)
But Σw_y(L_sm) = L_sm Σw_y(1). Eq(1) is solvable for L_sm
and
Σw_s(d_smy) = 0
or Σw_s(G_my) = Σw_y(T_smy - L_sm)......Eq(2)
But Σw_s(G_my) = G_my Σw_s(1). Eq(2) is solvable for G_my

The naive model is in fact a first step. Start with G=0. Solve (1) gives the normals, and solve (2) gives the resulting anomaly average.

It's not right, because as said, with non-zero G, (1) is not satisfied. But you can solve again, with the new G. Then solve (2) again with the new L. Repeat. When (and if) it has converged, both equations are satisfied.

Code

Here is some R code that I use. It starts with x and w arrays (dimensioned ns×nm×ny×, nm=12,nsm=ns*nm) as defined above:

G=0; wx=w*x; wr=rowSums(w)+1e-9; wc=colSums(w);
for(i in 1:8){
 L=rowSums(wx-w*rep(G,each=nsm))/wr
 G=colSums(wx-w*rep(L,ny))/wc
}

That's it. The rowSums etc do the averaging. I add 1e-9 to wr in case stations have no data for a given month. That would produce 0/0, which I want to be 0. I should mention that I have set all the NAs in x to 0. That's OK, because they have zero weight. Of course, in practice I monitor convergence and end the loop when achieved. We'll see that later.

In words

1. Initially, the offsets L are the long term station monthly averages
2. Anomalies are formed by subtracting L, and the global G by spatially averaging (carefully) the anomalies. Thus far, the naive method.
3. Form new offsets by subtracting the global anomaly from measured T, and averaging that over years.
4. Use the new offsets to make a spatially averaged global anomaly
5. Repeat 3 and 4 till converged.

A simple example

Here's another simple example showing that just taking anomalies over station data isn't enough. We'll forget months; the climate is warming 1C/year (change the units if you like), and there are two stations which measure this exactly over five years. But they have missing readings. So the x array looks like this:

Station A	1	2	3	-	-
Station B	-	-	3	4	5

and w looks like

w	1	1	1	0	0
	0	0	1	1	1

Here are the first four iterations. I've calculated the trend of G, which should be 1, and the weighted sum of squares of residuals.

	j
Iter	Trend	Sum squares
1	0.4	2
2	0.6	0.8889
3	0.7333	0.3951
4	0.8222	0.1756

The first row is the naive recon. The trend is way too low, and the sum squares is far from zero. Iteration shows gradual improvement. The columns do converge to 1 and 0 respectively.

The convergence speed depends on the overlap between stations. If there had not been overlap in the middle, the algorithm would have happily reported a zig-zag with the anomaly of each station as the global value. If I had extended the range of each station by 1 year, convergence is rapid. A real example would have lots of overlap. Here are plots of the comvergence:

A real example

Now I'll take x and w from a real TempLS run, from 1899 to April 2015. They are 9875×12×117 arrays, and returns 12*117 monthly global anomalies G. I'll give the complete code below. It takes four seconds on my PC to do six iterations:

Iter	Trend	Sum squares
1	0.5921	1
2	0.6836	0.994423
3	0.6996	0.994266
4	0.7026	0.99426
5	0.7032	0.99426
6	0.7033	0.99426

Again, the first row, corresponding to the naive method, significantly understates the whole-period trend. The sum of squares does not change much, because most of the monthly deviation is not explained by global change. The weights I used were mesh-based, and as a quick check with the regular results

	Nov	Dec	Jan	Feb	Mar	Apr
Regular run 18 May	0.578	0.661	0.662	0.697	0.721	0.63
This post	0.578	0.661	0.662	0.697	0.721	0.630

To 3 figures, it is the same.

More code

So here is the R code including preprocessing, calculation and diagnostics. I've used color to show these various parts:

for(i in 1:8){
load("TempLSxw.sav") # x are raw temperatures, w wts from TempLS
x[is.na(x)]=0
d=dim(x); # 9875   12  117
L=v=0;  
wr=rowSums(w,dims=2)+1e-9; # right hand side
wc=colSums(w)+1e-9; # right hand side
for(i in 1:6){
 L=rowSums(w*x-w*rep(G,each=d[1]),dims=2)/wr  # Solve for L
 G=colSums(w*x-w*rep(L,d[3]))/wc # then G
 
 r=x-rep(L,ny)-rep(G,each=ns)  # the residuals
 v=rbind(v,round(c(i,
 lm(c(G)~I(1:length(G)))$coef[2]*1200,sum(r*w*r)),4))
}

colnames(v)=c("Iter","Trend","Sum sq re w")
v[,3]=round(v[,3]/v[2,3],6)  # normalising SS
G=G-rowMeans(G[,62:91+1]) # to 1961-90 anomaly base
write(html.table(v),file="tmp.txt")

Emissions hiatus?

John Baez has a post on the latest preliminary data from EIA for 2014. Global CO2 emissions are the same as for 2013 at 32.3 Gtons. EIA says that it is the first non-increase for 40 years that was not tied to an economic downturn.

They attribute the pause to greater use of renewables, mentioning China. Greenpeace expands on this, saying that China's use of coal dropped by 8%, with a consequent 5% drop in CO2 emission. They give the calc with sources here. This source says April coal mined in China was down 7.4% on last year, which they do partly attribute to economic slowdown there.

It's just one year, and may be influenced by China's economy. We'll see.

JAXA Hiatus

If you're following JAXA Arctic sea ice, they have announced an interruption for maintenance, until May 20. Pity, it is getting interesting, with 2015 now at record low levels for the day. Maybe it was going off the rails - we'll see.

NSIDC NH isn't looking so reliable either - they showed a massive re-freeze yesterday.

Update - JAXA has unhelpfully taken down all its data, and replaced it with the warning (in Japanese), which overwrote my local data. So no JAXA table data here till that is restored. Fortunately the plot is OK.
Update - Jaxa is back.