Tamino over at  Open Mind has a new post detailing his approach for calculating temperature averages. See his post here.

His method is based on the Berkeley method as he notes and he uses it primarily for calculating regional or local temperature averages. Read his post for the math details behind the approach. I got my hands on the code last night and set about trying to figure out how to integrate his approach in the RghcnV3 package.  That, of course, meant I stopped work on version 1.2 which is ready to hit CRAN.  Tamino, posted his code in a pdf. That caused a little trouble with copying and pasting, but in the end it was easy to clean up.

His code would be a good example for other folks ( who shall go unmentioned) to follow. It was clean and understandable. I might complain a bit about variable names and ask for less terse names, but  since I got the general idea of the algorithm it was easy to follow.

Here is his code

# Here's the code:
##################################
##################################
# align
# combine station data
# t is time
# x is data
# part is station ID
# offres is resolution for offsets
# tround is rounding digits for t
# xround is rounding digits for x
##################################
##################################
align = function(t,x,part,full=F,tround=NULL,xround=4,offres=.001){
##########################
# make "part" a factor
# and use 1st as reference
##########################
part = as.factor(part)
nparts = length(levels(part))
if (nparts < 2){
zout = data.frame(t,x)
zout = zout[order(t),]
zout = list(data=zout,offsets=0)
} else{
#######################
# form giant data frame
#######################
for (jpart in 1:nparts){
tt = t[part == levels(part)[jpart]]
xx = x[part == levels(part)[jpart]]
xx = xx[order(tt)]
tt = tt[order(tt)]
if (length(tround) > 0){
tt = round(tt,tround)
}
z1 = data.frame(tt,xx)
names(z1) = c("t",levels(part)[jpart])
if (jpart == 1){
   zz = z1
} else{
    zz = merge(zz,z1,by=1,all=T)
}
}
###############################
# iteratively determine offsets
# and merged values
###############################
zdat = zz[,2:(nparts+1)]
zdat = as.matrix(zdat)
ntimes = dim(zz)[1]
mu = rep(0,nparts)
gam = rep(NA,ntimes)
dmu = 9999
while (dmu > offres){
for (jj in 1:ntimes){
xx = zdat[jj,]
xx = xx-mu
gam[jj] = mean(xx,na.rm=T)
}
oldmu = mu
for (jj in 1:nparts){
xx = zdat[,jj]
xx = xx-gam
mu[jj] = mean(xx,na.rm=T)
}
##################################
# shift offsets so 1st offset is 0
##################################
mu = mu-mu[1]
dmu = sum(abs(mu-oldmu))
}
############
# final mean
############
num = rep(NA,ntimes)
se = rep(NA,ntimes)
std.dev = rep(NA,ntimes)
for (jj in 1:ntimes){
xx = zdat[jj,]
xx = xx-mu
xx = xx[is.finite(xx)]
num[jj] = length(xx)
gam[jj] = mean(xx)
std.dev[jj] = sd(xx)
}
se = std.dev/sqrt(num)
zout = data.frame(x=gam,num,se,std.dev)
if (length(xround) > 0){
zout = round(zout,xround)
}
zout = data.frame(t=zz[,1],zout)
if (full){
zd = data.frame(zdat)
names(zd) = levels(part)
zout = data.frame(zout,zd)
}
zout = list(data=zout,offsets=mu)
}
zout
}