Using RghcnV3 - a very simple TempLS

I've been experimenting with Steven Mosher's R package RghcnV3. It has many useful features, and I think the R package format serves it well. I'm not as enthusiastic as Steven about the zoo and raster structures, at least for this purpose, but still, I may yet be convinced.

Anyway, it certainly gets some of the messy initial data structuring out of the way. So since I am currently working on Ver 2.2, I thought I would put together a very simple version using some of the more recent ideas, in the RghcnV3 style.
Update - I've added a version of the code with detailed comments at the end.

TempLS basics

Firstly, a reminder of how TempLS works. After data rearrangement, we have an array of temperature readings, over some number S of stations, and over M months (usually 12, but can do seasons etc) and Y years. This is a SxMxY array with many missing values.

I'll explain the math in a similar style to this account of V2. I'll use a notation, described there, in which arrays (matrices etc) are represented with suffixes, and equations including them are true for the natural range of values. To save summation signs, I use a summation convention in which a repeated index is assumed to be summed over the range. But to emphasise that, I use colors. Summed indices are blue, independent ones are red. There are some I need to exempt from the summation rule, and these are green.

So we assume that the observed temps x can be modelled with

• a local temperature L, which varies over station and months, but not over years. Think of it as a monthly average, and
• a global temperature G uniform over months and stations, but varying with years. This corresponds to the global anomaly that we are after.

So we write a model equation

x_smy  ~  L_sm  +  G_y
I'll repeat here the least squares part from that earlier post. It's useful to maintain a convention that every term has to include every red index. So I'll introduce an identity I_y, for which every component is 1, and rewrite the model as
x_smy  ~  I_yL_sm  +  I_smG_y

Weighted least squares.

The L and G are parameters to be fitted. We minimise a weighted least squares expression for the residuals. The weight function w has two important roles.

• It is zero for missing values
• It needs to ensure that stations make "fair" contributions to the sum, so that regions with many stations aren't over-represented. This means that the sum with weights should be like a numericam integral over the surface. w should be inversely proportional to density.

So we want to minimise w_(smy)(x_smy - I_yL_sm - I_smG_y)²

Differentiation

This is done by differentiating this expression wrt each parameter component. In the summation convention, when you differentiate a paired index, it disappears. The remaining index goes from blue to red. The original scalar SS becomes a system of equations, as the red indices indicate:

Minimise:	w(smy)(xsmy - IyLsm - IsmGy)2
Differentiate wrt L:	w(smy)Iy(xsmy - IyLsm - IsmGy)
Differentiate wrt G:	w(smy)Ism(xsmy - IyLsm - IsmGy)

Notice what we have here - linear equations in L and G, and all we need is the data x and the weights w.

Equations to be solved

These are the equations to be solved for L and G. It can be seen as a block system:
A*L + B*G  = X1
B^T*L + C*G  = X2
combining s and m as one index, and where ^T signifies transpose. Examining the various terms, one sees that A and C are diagonal matrices, so the overall system is symmetric.

The first equation can be used to eliminate L, so
(C -  B^TA^-1B)*G = X2  -  B^TA^-1X1
This is not a large system - G has of order 100 components (years).

Implementation with RghcnV3

There are just three steps. We have to get the array x. RghcnV3 can get us (as explained in the manual as far as the output of readV3Data(). This is a large array ordered first by months (in rows), then by year, then by station. There is then a function v3ToZoo() in V1.3) which should get it into a structured array. Unfortunately, I ran out of memory at that stage. But the further rearrangement was a basic step in TempLS, so I used the sortdata() routine to do that. The arguments are the data array, and a vector with the start and end year of the period you want to look at.The resulting x array is ordered as suggested above, by station, month and year.

Weighting

Then we need a weight function w. A basic method is to divide the surface into 5x5 lat/lon cells, and get a station density by dividing the area of each cell by the number its stations reporting in each month. This is done in the routine makeweights, which returns the array w.

Solving

So that's it. We have the needed w and x; the rest is just some summing over rows etc to get the block matrices, and then a block inversion. This is done in the routine solveTemperature, which returns the temperature anomalies about a mean zero. I've added a plot statement.

The code

As promised, just three routines. Start with dat, the array from readV3Data(), and your choice of years yrs[1:2]. There is another version with detailed comments below:

sortdata = function(v,yrs=c(1900,2010)) # Sorts data from v2.mean etc into matrix array x
{
nd=dim(v)[1];
y0=yrs[1]-1;
iq=c(0,which(diff(v[,1])>0),nd);
jq=length(iq)-1;
x=array(NA,c(jq,12,yrs[2]-y0));
for(i in 1:jq){ # loop over stations
kq=c(iq[i]+1,iq[i+1]);
iy=as.matrix(v[kq[1]:kq[2],]);
oy=(iy[,2] > y0) & (iy[,2]<=yrs[2]);
if(sum(oy)<2)next;
iz= iy[oy,];
x[i,,iz[,2]-y0] = t(iz[,3:14]);
}
x;
}

Then the routine to make the weights - needs x and lat/lon:
makeweights = function( x, inv)
{
d = dim(x); wt = array(NA,d);
ns = floor(inv$Lat/5)*72 + floor(inv$Lon/5) + 1333;
wtt = sin((1:36-0.5)*pi/36);
w_area = rep(wtt,each=72);
for(i in 1:d[2]){
for(j in 1:d[3]){
ow = !is.na(x[,i,j]);
n=ns[ow];
kc=tabulate(c(n,72*36));
ok=kc>0; wc=w_area;
wc[ok]=w_area[ok]/kc[ok];
wt[ow,i,j]= wc[n];
}
}
wt;
}

Finally, the solve routine:

solvetemperature=function(wt,x)
{
wt[is.na(wt)]=0;x[is.na(x)]=0;
x=x*wt; d=dim(x); d=c(d[1]*d[2],d[3]);
dim(x)=dim(wt)=d;
A=rowSums(wt)+1.0e-9;
R1=rowSums(x);
C=colSums(wt);
R2=colSums(x);
aw=array(rep(1/A,d[2])*wt,d);
b = drop(R2-R1 %*% aw);
S=diag(C) - t(wt) %*% aw;
S[1,1]=S[1,1]+1; # there is a dof to be fixed
y=solve(S,b);
y-mean(y);
}

And the program to run it all:

x=sortdata(dat,yrs);
wt=makeweights(x,inv);
T=solvetemperature(wt,x);
png("TempLSv0.png",width=500);
plot(yrs[1]:yrs[2],T,type="l");
dev.off()

And a (no-frills) result:





Here is a version of the code above with line by line comments.

#You'll want your own version of yrs
sortdata = function(v,yrs=c(1900,2010)) # Sorts data from v2.mean etc into matrix array x.
#v is just the v3.mean file read in - one year per row, one station block after another. Within a block, the years are in order.
{

   nd=dim(v)[1]; # Number of stations in inventory
   y0=yrs[1]-1; # Year zero
   iq=c(0,which(diff(v[,1])>0),nd); # Endpoints of each station block of data (where the stat num ber increments)
   jq=length(iq)-1; # Number of stations in data set
   x=array(NA,c(jq,12,yrs[2]-y0)); # setting up the big temp matrix
   for(i in 1:jq){ # loop over stations
       kq=c(iq[i]+1,iq[i+1]); # marks the block of station i
       iy=as.matrix(v[kq[1]:kq[2],]); # copies that block of info
       oy=(iy[,2] > y0) & (iy[,2]<=yrs[2]); # selects the years within range
       if(sum(oy)<2)next; # reject stations with less than a year of data (can't get an anomaly with 1 dof)
       iz= iy[oy,]; # the selected block
       x[i,,iz[,2]-y0] = t(iz[,3:14]); # copy into x
   }
x; # return x
}#;
makeweights = function( x, inv)
# To return wt, a matrix with the same structure as x
#ie a weight for every station/month/year, but zero or NA if missing
{
  d = dim(x); wt = array(NA,d); # set up the big blank
  ns = floor(inv$Lat/5)*72 + floor(inv$Lon/5) + 1333;
# ns lists a (5x5 degree) cell number for every station in the inventory
  wtt = sin((1:36-0.5)*pi/36); w_area = rep(wtt,each=72);
# w_area is an area weight for each cell
  for(i in 1:d[2]){ # loop over months
    for(j in 1:d[3]){ # loop over years
       ow = !is.na(x[,i,j]); # which stations reported in i/j
       n=ns[ow]; # and what were their cell numbers
       kc=tabulate(c(n,72*36)); # frequency of cell numbers
       ok=kc>0; # logical for empty cells (that month)
       wc=w_area; #Initialize the weights to the areas
       wc[ok]=w_area[ok]/kc[ok]; # now divide the weights by number of stats, to get inverse density
       wt[ow,i,j]= wc[n]; # set the weights for that i,j
    }
   }
 wt; # return the big weight matrix
}
#Now we have the weights wt and temp matrix x, can solve linear systen for temp anomaly
solvetemperature=function(wt,x)
{
# to avoid NA fuss, missing readings have zero weight, so temp doesn't matter
   wt[is.na(wt)]=0;x[is.na(x)]=0;
   # redimension to a matrix with s&m vs year
   x=x*wt; d=dim(x); d=c(d[1]*d[2],d[3]);
   dim(x)=dim(wt)=d;
   A=rowSums(wt)+1.0e-9; # Top left diag block matrix
# There can be missing lines with all zero. Instead of eliminating them, just add something to the diag to avoid a 0/0 error. No effect on result.
   R1=rowSums(x); # Corresponding RHS
   C=colSums(wt); #Bottom right diag matrix
   R2=colSums(x); # Corresponding RHS
   aw=array(rep(1/A,d[2])*wt,d); # LHS term after dividing by A to eliminate L
   b = drop(R2-R1 %*% aw); # Elimination op for RHS
   S=diag(C) - t(wt) %*% aw; # Local temp has now been eliminated, along with variables s and m.
#Now the equation is just in yearly T anomalies over years. Symmetric.
# S is singular, because an arbitrary constant could be added to G and subtracted from L.
# So arbitrarily set anomaly 1 to zero
   S[1,1]=S[1,1]+1;
# Now solve - S could be up to 150x150 say (yr[2]-yr[1])
   y=solve(S,b);
# Now reset to mean zero (fixing that arbitrariness)
   y-mean(y); # zero mean anomaly returned.
}
# End functions.
# Execution You start with v from data file (readV3Data()) and inv from inventory (readInventory())
# All you need from inventory is lat/lon
x=sortdata(dat,yrs);
# Many other makeweights() are possible
wt=makeweights(x,inv);
T=solvetemperature(wt,x);
# A simple plot on a png file; comment out to see on screen (also dev.off())
png("TempLSv0.r");
plot(yrs[1]:yrs[2],T,type="l",xlab="Year");
dev.off();