Ryan’s code for RegEM: The Magic Teleconnection Machine

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

### SEGMENT: PLOTTING FUNCTIONS ###

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

###

### 1. set.colors: Generates a color map with blues for negative values and

### reds for positive values. Also can set a specified range to a single

### color for highlighting. Returns a vector of colors, starting with

### -limit and ending with +limit.

###

### rng: The range between 0 and +limit or -limit. Range is ALWAYS

### a positive number (negative numbers will invert the color scale

### without inverting the labels).

###

### divs: The number of different colors between (and not including)

### zero and +limit (and, by extension, between -limit and zero).

###

### hlite: The single value (hlite=number) or range (hlite=c(number1, number2)

### of values to be plotted with color “hcol”.

###

### hcol: The color used for highlighting points within the specified range.

###

### sides: 0: Scale is -rng to +rng

### 1: Scale is 0 to +rng

### -1: Scale is -rng to 0

###

### 2. set.legend: Generates the automatic legend for plots utilizing a color

### map.

###

### x, y: The x and y coordinates for the upper left corner of the legend box.

###

### sl, st, spl: Scaling factors for the box, legend text, and size of the

### color scale bar, respectively.

###

### h: The highlighted range, used for generating the x-axis label describing

### the range.

###

### rng: The range between 0 and +limit or -limit. Range is ALWAYS

### a positive number (negative numbers will invert the color scale

### without inverting the labels).

###

### ml: A list containing the main title and subtitle for the plot.

###

### mbx, mxs: Logical values for plotting a box around the plot and the plot

### x-y axes, respectively.

###

### mc: Vector containing the colors for the color scale.

###

### divs: The number of colors between (but not including) zero and +limit.

###

### s.main, s.sub, s.axis: Scaling factors for the main title, subtitle, and

### axis labels, respectively.

###

### sides: 0: Scale is -rng to +rng

### 1: Scale is 0 to +rng

### -1: Scale is -rng to 0

###

### 3. plt.map(dat, coord, labels, cols, legends, norm, sp, …): Plots points

### on a projection of Antarctica.

###

### dat: The data to be plotted.

###

### coord: Lon/Lat coordinates of the points. Longitude should be in column

### 1 and latitude in column 2.

###

### cols: Vector of colors. Default is NA, which results in automatic

### generation of a color scale.

###

### legends: Default is NA, which results in automatic generation of a legend

### for the color scale. If manual legends are desired, legends

### requires a list of:

### [[1]] Legend location (i.e., “bottomleft”)

### [[2]] Vector of legend text

### [[3]] Value for pch

### [[4]] Vector of colors corresponding to legend text

### [[5]] cex value for legend text and points

### [[6]] Background color

###

### norm: If TRUE, normalizes the plotted data to a scale of +/- 1.

###

### sp: Size parameter for plotted points.

###

### … : Parameters to be passed to set.colors and set.legend for automatic

### color map and legend generation.

###

set.colors=function(rng, divs, hlite, hcol, sides) {

### Make color scale

cold.colors=matrix(nrow=divs); hot.colors=matrix(nrow=divs)

for(i in 1:divs) {

lum.ratio=25*(i/divs)

c1=ifelse(lum.ratio<15, 255-(17*lum.ratio), 0)

c2=ifelse(lum.ratio<15, 255, 255-(10*(lum.ratio-15)))

h1=ifelse(255-(12*lum.ratio)>0, 255-(12*lum.ratio), 0)

h2=ifelse(255-(25*lum.ratio)>0, 255-(25*lum.ratio), 0)

cold.colors[divs+1-i]=rgb(red=c1, green=c1, blue=c2, maxColorValue=255)

hot.colors[i]=rgb(red=c2, green=h1, blue=h2, maxColorValue=255)

}

### Define the zero point color. If few enough divs for a space between each color,

### then use gray for the zero point color. Otherwise, use white.

neut.col=ifelse(divs<11, “#EEEEEE”, “#FFFFFF”)

### Combine the colors into the color scale vector

if(sides==0) {

mapcolors=c(“darkmagenta”, as.vector(cold.colors), neut.col, as.vector(hot.colors), “black”) }

if(sides==-1) {

mapcolors=c(“darkmagenta”, as.vector(cold.colors), neut.col) }

if(sides==1) {

mapcolors=c(neut.col, as.vector(hot.colors), “black”) }

### If hlite is used, replace appropriate portions of the color scale with hcol

if(!is.na(hlite[1])) {

hlite[2]=ifelse(is.na(hlite[2]), hlite[1], hlite[2])

start.pos=ifelse(sides==0, (divs+2), 1)

start.pos=ifelse(sides==-1, length(mapcolors), start.pos)

highlmin=start.pos+(round((divs+1)*(min(hlite[1], hlite[2])/(rng))))

highlmin=ifelse(highlmin<1, 1, highlmin)

highlmax=start.pos+(round((divs+1)*(max(hlite[1], hlite[2])/(rng))))

highlmax=ifelse(highlmax>length(mapcolors), length(mapcolors), highlmax)

mapcolors[highlmin:highlmax]=hcol

}

### Return the color scale vector

mapcolors

}

set.legend=function(x, y, sl, st, spl, h, hcol, rng, ml, mbx, mxs, mc, divs, s.main, s.sub, s.axis, sides) {

### Draw the box for the legend

rect(xleft=x, xright=x+(.6*sl), ytop=y, ybottom=y-(.080*sl), col=”white”)

### If sides=0, plot the neutral color (exact center of the color scale)

if(sides==0) {points(x=x+(.3*sl), y=y-(.023*sl), col=mc[2+divs], pch=15, cex=1.3*spl)}

### Plot the rest of the color scale

if(sides==-1 | sides==1) {

for(i in 1:divs) {

points(x=x+.1*sl+(.4/(divs+1))*i*sl, y=y-(.023*sl), col=mc[1+i], pch=15, cex=1.3*spl)

}

}

if(sides==0) {

for(i in 1:divs) {

points(x=x+(.3*sl)-(.2/(divs+1))*i*sl, y=y-(.023*sl), col=mc[divs+2-i], pch=15, cex=1.3*spl)

points(x=x+(.3*sl)+(.2/(divs+1))*i*sl, y=y-(.023*sl), col=mc[divs+2+i], pch=15, cex=1.3*spl)

}

}

### Plot the -limit and +limit colors of the color scale

points(x=x+(.1*sl), y=y-(.023*sl), col=mc[1], pch=15, cex=1.3*spl)

points(x=x+(.5*sl), y=y-(.023*sl), col=mc[length(mc)], pch=15, cex=1.3*spl)

### Plot the units, -limit label, and +limit label

if(sides==0) {lmt1=-rng;lmt2=rng}

if(sides==-1) {lmt1=-rng; lmt2=0}

if(sides==1) {lmt1=0;lmt2=rng}

text(x=x+(.3*sl), y=y-(.055*sl), ml[1], cex=st)

text(x=x+(.085*sl), y=y-(.055*sl), lmt1, cex=st)

text(x=x+(.51*sl), y=y-(.055*sl), lmt2, cex=st)

### If hlite is used, define the text describing the highlighted range

if(!is.na(h[2])) {

ltext=paste(min(h[1],h[2]), “to”, max(h[1], h[2]))

ltext=ifelse(min(h[1],h[2])<=lmt1, paste(max(h[1],h[2]), “and Less”), ltext)

ltext=ifelse(max(h[1],h[2])>=lmt2, paste(min(h[1],h[2]), “and Greater”), ltext)

}

rtext=ifelse(!is.na(h[2]), ltext, h[1])

xtext=ifelse(!is.na(h[1]), paste(rtext, ml[1], “Highlighted”), NA)

### Print the title, subtitle, and highlighted range text

title(main=ml[2], sub=ml[3], xlab=xtext, cex.main=s.main, cex.sub=s.sub, cex.axis=s.axis)

### Draw box and/or axes?

if(mbx==TRUE) {box()}

if(mxs==TRUE) {map.axes()}

}

map.lbls=c(“Deg C/Decade”, “”, “”)

plt.map=function(dat, coord=sat.coord, rng=.5, divs=1000, labels=map.lbls, norm=FALSE, hlite=NA, hcol=”#99FF33″, mbx=TRUE, mxs=FALSE, x=-.3, y=-.35, sl=1, st=1, sp=1, spl=1, s.main=1, s.sub=1, s.axis=1, cols=NA, legends=NA, sides=0) {

library(mapproj)

### Set up Antarctic projection and lat/lon lines

map(“world”, proj=”azequidistant”, orient=c(-90, 0, 0), ylim=c(-90, -60))

a=seq(-180, 150, 30)

for(i in 1:length(a)) {lines(mapproject(list(x=rep(a[i], 181), y=c(-90:90))), col=4, lty=2)}

a=c(-180,0)

for(i in 1:length(a)) {lines(mapproject(list(x=rep(a[i], 181), y=c(-90:90))), col=4, lty=1)}

a=seq(-50,-80,-10)

for(i in 1:4) {lines(mapproject(list(y=rep(a[i], 72), x=seq(-177.5, 177.5,5))), col=4, lty=2)}

### If no user provided color scale, make the color scale

if(is.list(cols)) {mapcolors=cols} else {mapcolors=set.colors(rng, divs, hlite, hcol, sides)}

### Normalize the data to a maximum of +/- 1?

nm=ifelse(norm==TRUE, max(abs(as.vector(dat))), 1)

### Assign the color scale

if(!is.list(cols)) {

start.pos=ifelse(sides==0, (divs+2), 1)

start.pos=ifelse(sides==-1, length(mapcolors), start.pos)

temp=start.pos+(round((divs+1)*(dat/(rng*nm))))

temp=ifelse(temp<1, 1, temp)

temp=ifelse(temp>length(mapcolors), mapcolors[length(mapcolors)], mapcolors[temp])

} else {temp=cols}

### Plot points

xy=mapproject(coord[, 1], coord[, 2])

points(xy$x, xy$y, pch=15, cex=sp, col=temp)

### Overlay map

map(“world”, proj=”azequidistant”, orient=c(-90, 0, 0), ylim=c(-90, -60), fill=FALSE, add=TRUE)

### Set up legend and labels

if(!is.list(legends)) {

legends=set.legend(x, y, sl, st, spl, h=hlite, hcol, rng, ml=labels, mbx, mxs, mc=mapcolors, divs, s.main, s.sub, s.axis, sides)

} else {

legend(legends[[1]], legends[[2]], pch=legends[[3]], col=legends[[4]], cex=legends[[5]], bg=legends[[6]])

title(main=labels[1], sub=labels[2])

}

}

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

calc.anom=function(dat) {

### A present from Steve

for(i in 1:12) { sequ = seq(i, nrow(dat), 12)

dat[sequ, ] = scale(dat[sequ, ], scale=F) }

dat

}

ssq=function(x) {sum(x^2)}

regem.pca=function(X, maxiter=999, tol=.01, regpar=3, method=”default”, split=FALSE, l=0) {

### Define regem, which contains the output of each iteration

regem=list()

### Initial infill of NAs with zeros

n=nrow(X);p=ncol(X)

indmis=is.na(X)

nmis=sum(indmis)

dofC=n-1

trunc=min(regpar, n-1, p)

dofS=dofC-trunc

### Center and infill NAs with zeros

X=scale(X, scale=FALSE)

M=attributes(X)$”scaled:center”

X[indmis]=0

### Define loop start

iter=0

dXmis=10^6

while(iter<=maxiter & dXmis>tol) {

iter=iter+1

peff_ave=0

### Define cov(X)

C=t(X) %*% X/dofC

CovRes=array(0, dim=c(p, p))

### Define sd

D=sqrt(diag(C))

### Normalize to sd=1

X=X %*% diag(1/D)

C=diag(1/D )%*% C %*% diag(1/D)

### Find the singular value decomposition of the input matrix

### using regpar to truncate the number of right/left singular

### vectors to be computed and calculate the truncated rank matrix.

svd.X= svd(X, nu=regpar, nv=regpar)

if(regpar>1) {

Xmis= svd.X$u %*% diag(svd.X$d[1:regpar]) %*% t(svd.X$v)

}

if(regpar==1) {

Xmis= as.matrix(svd.X$u) %*% as.matrix(svd.X$d[1]) %*% t(svd.X$v)

}

### Rescale to original scaling

X=X1a=scale(X, center=FALSE, scale=1/D)

Xmis=scale(Xmis, center=FALSE, scale=1/D)

C= diag(D) %*% C %*% diag(D)

CovRes=diag(D) %*% CovRes %*% diag(D)

dXmis=sqrt(ssq(Xmis[indmis]-X[indmis]))/sqrt(nmis)

### Update and recenter

X[indmis]=Xmis[indmis]

X=scale(X, scale=FALSE)

Mup=attributes(X)$”scaled:center”

M=M+Mup

C=t(X) %*% X/dofC

### Store the iteration results and end the iteration

if(method==”verbose”) {

regem[[iter]]=list(X=X, C=C, B=B, S=S, dXmis=dXmis) } else {

regem[[iter]]=list(X=X, dXmis=dXmis)

}

}

print(“”, quote=FALSE); print(“RegEM Performance”, quote=FALSE); print(“=================”, quote=FALSE)

if(split==TRUE) {print(paste(“PC”, l, “Statistics:”), quote=FALSE);print(“—————-“, quote=FALSE)}

### Inform if results did not converge

if(dXmis>tol) {print(paste((iter-1), “iterations reached without convergence: dXmis =”, dXmis), quote=FALSE)}

### Inform number of iterations if results converged

if(dXmis<=tol) {print(paste((iter), “iterations for convergence: dXmis =”, dXmis), quote=FALSE)}

### Return the list of iteration results

N=length(regem)

regem[[N]]$svd=svd.X

regem[[N]]$Xmis=Xmis

regem

}

get.PCs=function(dat, base, PC.max=3, st=1982) {

dat=t(dat)

dat.svd=svd(dat)

### Multiply the right singular vectors through and store as a time series

PCs=ts(t((dat.svd$d)*t(dat.svd$v)), start=st, freq=12)

### Print explained variance

var.tot=0

print(“PC Selection”, quote=FALSE); print(“============”, quote=FALSE)

print(“Explained variance:”, quote=FALSE); print(“”, quote=FALSE)

for(i in 1:PC.max) {

print(paste(“PC “, i, “: “, round((dat.svd$d[i]^2/ssq(dat.svd$d))*100,2), “%”, sep=””), quote=FALSE)

var.tot=var.tot+(dat.svd$d[i]^2/ssq(dat.svd$d))*100

}

print(“”, quote=FALSE); print(paste(“Total explained variance: “, round(var.tot,2), “%”, sep=””), quote=FALSE); print(“”, quote=FALSE)

### Return the desired number of PCs

PCs=list(ts.union(base, PCs[, 1:PC.max]), dat.svd$u[, 1:PC.max], dat.svd$v, dat.svd$d)

PCs

}

get.recon=function(dat.predict, dat.cov, PC.max) {

if(PC.max==1) {

recon=array(0, dim=c(600, length(dat.cov)))

for(i in 1:length(dat.cov)){

recon[, i]=ts(dat.predict*dat.cov[i], start=1957, freq=12)

} } else {

recon=array(0, dim=c(600, nrow(dat.cov)))

for(i in 1:nrow(dat.cov)) {

temp.rec=ts(t(t(dat.predict)*(dat.cov[i, ])), start=1957, freq=12)

recon[,i]=ts(rowSums(temp.rec), start=1957, freq=12)

}

}

recon

}

stat.trends=function(dat, st, en) {

dat=window(dat, st, en)

### Store trends as deg C per decade

temp=lm(dat~I(time(dat)))[[1]][2, ]*10

temp

}

info=c(“”, “Elements in returned list:”,”=========================”, “[[1]]: Reconstruction time series”,”[[2]]: Reconstruction overall monthly means time series”,”[[3]]: Imputed PCs”,”[[4]]: Final RegEM iteration (if split RegEM, a list containing all of the final iterations)”,”[[5]][[1]]: Ground station and PCs matrix (input to RegEM)”,”[[5]][[2]]: Left singular vectors (1 – PC.max)”,”[[5]][[3]]: Right singular vectors (all)”,”[[5]][[4]]: Eigenvalues (all)”, “[[6]]: Linear trends (1957-2006)”, “[[7]]: Linear trends (1965-2006)”, “[[8]]: Linear trends (1982-2006)”, “[[9]]: Linear trends (1982-1994.5)”, “[[10]]: Linear trends (1994.5-2006)”)

do.recon=function(dat=avhrr.anom, base=reader.steig.anom, PC.max=3, maxiter=999, tol=.01, regpar=41, method=”default”, base.st=1957, dat.st=1982, rec.en=c(2006,12), infor=info, tr=TRUE) {

PCs=get.PCs(dat, base, PC.max, dat.st)

PC.dat=PCs[[1]]

PC.cov=PCs[[2]]

regem=regem.ftls(PC.dat, maxiter, tol, regpar, method)

regem.last=regem[[length(regem)]]$X

regem.PCs=regem.last[, -(1:ncol(base))]

recon=ts(get.recon(regem.PCs, PC.cov, PC.max), start=base.st, freq=12)

trend.all=0; trend.65=0; trend.sat=0; trend.82=0; trend.94=0

if(tr==TRUE) {

trend.all=stat.trends(recon, 1957, c(2006,12))

trend.65=stat.trends(recon, 1965, c(2006,12))

trend.sat=stat.trends(recon, 1982, c(2006,12))

trend.82=stat.trends(recon, 1982, c(1994,6))

trend.94=stat.trends(recon, c(1994,7), c(2006,12))

}

recon.means=ts(rowMeans(recon), start=base.st, freq=12)

all.recon=list(recon, recon.means, regem.PCs, regem.last, PCs, trend.all, trend.65, trend.sat, trend.82, trend.94)

print.table(info)

all.recon

}

do.split.recon=function(dat=avhrr.anom, base=reader.steig.anom, PC.max=3, maxiter=999, tol=.01, regpar=3, method=”default”, base.st=1957, dat.st=1982, rec.en=c(2006,12), infor=info, tr=TRUE) {

PCs=get.PCs(dat, base, PC.max, dat.st)

PC.cov=PCs[[2]]

regem.last=list()

regem.PCs=matrix(nrow=nrow(base), ncol=PC.max)

for(i in 1:PC.max) {

PC.dat=PCs[[1]][, c(1:ncol(base), ncol(base)+i)]

regem=regem.pca(PC.dat, maxiter, tol, regpar, method, split=TRUE, l=i)

regem.last[[i]]=regem[[length(regem)]]$X

regem.PCs[, i]=regem.last[[i]][, -(1:ncol(base))]

}

recon=ts(get.recon(regem.PCs, PC.cov, PC.max), start=base.st, freq=12)

trend.all=0; trend.65=0; trend.sat=0; trend.82=0; trend.94=0

if(tr==TRUE) {

trend.all=stat.trends(recon, 1957, c(2006,12))

trend.65=stat.trends(recon, 1965, c(2006,12))

trend.sat=stat.trends(recon, 1982, c(2006,12))

trend.82=stat.trends(recon, 1982, c(1994,6))

trend.94=stat.trends(recon, c(1994,7), c(2006,12))

}

recon.means=ts(rowMeans(recon), start=base.st, freq=12)

all.recon=list(recon, recon.means, regem.PCs, regem.last, PCs, trend.all, trend.65, trend.sat, trend.82, trend.94)

print.table(info)

all.recon

}

### Load files

lon=read.table(“tir_lons.txt”)

lat=read.table(“tir_lats.txt”)

sat.coord = cbind(lon,lat)

load(“reader.data.04-04-2009.RData”)

load(“reader.steig.04-12-2009.RData”)

grid=scan(“temp.dat”, n= -1)

avhrr=ts(t(array(grid, dim=c(5509, 300))), start=1982, freq=12); rm(grid)

grid=scan(“ant_recon.txt”, n= -1)

ant.recon=ts(t(array(grid, dim=c(5509, 600))), start=1957, freq=12); rm(grid)

### Calculate anomalies

reader.anom=calc.anom(reader.data); reader.steig.anom=calc.anom(reader.steig); avhrr.anom=calc.anom(avhrr)

### Get AVHRR and main TIR recon trends

# Use to calculate:

# m.trends=list()

# a.trends=list()

# m.trends[[6]]=stat.trends(ant.recon, 1957, c(2006,12))

# m.trends[[7]]=stat.trends(ant.recon, 1965, c(2006,12))

# m.trends[[8]]=stat.trends(ant.recon, 1982, c(2006,12))

# m.trends[[9]]=stat.trends(ant.recon, 1982, c(1994,6))

# m.trends[[10]]=stat.trends(ant.recon, c(1994,7), c(2006,12))

# a.trends[[6]]=0

# a.trends[[7]]=0

# a.trends[[8]]=stat.trends(avhrr.anom, 1982, c(2006,12))

# a.trends[[9]]=stat.trends(avhrr.anom, 1982, c(1994,6))

# a.trends[[10]]=stat.trends(avhrr.anom, c(1994,7), c(2006,12))

# Alternately load saved trends:

load(“avhrr.trends.RData”)

load(“TIR.trends.RData”)

m.means=ts(rowMeans(ant.recon), start=1957, freq=12)

a.means=ts(rowMeans(avhrr.anom), start=1982, freq=12)

### Do a 3-PC, 3-regpar recon:

#recon=do.recon()

=============================================================

=============================================================####################################
### SEGMENT:  PLOTTING FUNCTIONS ###
####################################
###
###  1.  set.colors:  Generates a color map with blues for negative values and
###      reds for positive values.  Also can set a specified range to a single
###      color for highlighting.  Returns a vector of colors, starting with
###      -limit and ending with +limit.
###
###      rng:  The range between 0 and +limit or -limit.  Range is ALWAYS
###            a positive number (negative numbers will invert the color scale
###            without inverting the labels).
###
###      divs:  The number of different colors between (and not including)
###             zero and +limit (and, by extension, between -limit and zero).
###
###      hlite:  The single value (hlite=number) or range (hlite=c(number1, number2)
###              of values to be plotted with color “hcol”.
###
###      hcol:  The color used for highlighting points within the specified range.
###
###      sides:  0:  Scale is -rng to +rng
###              1:  Scale is 0 to +rng
###             -1:  Scale is -rng to 0
###
###  2.  set.legend:  Generates the automatic legend for plots utilizing a color
###      map.
###
###      x, y:  The x and y coordinates for the upper left corner of the legend box.
###
###      sl, st, spl:  Scaling factors for the box, legend text, and size of the
###                    color scale bar, respectively.
###
###      h:  The highlighted range, used for generating the x-axis label describing
###          the range.
###
###      rng:  The range between 0 and +limit or -limit.  Range is ALWAYS
###            a positive number (negative numbers will invert the color scale
###            without inverting the labels).
###
###      ml:  A list containing the main title and subtitle for the plot.
###
###      mbx, mxs:  Logical values for plotting a box around the plot and the plot
###                 x-y axes, respectively.
###
###      mc:  Vector containing the colors for the color scale.
###
###      divs:  The number of colors between (but not including) zero and +limit.
###
###      s.main, s.sub, s.axis:  Scaling factors for the main title, subtitle, and
###                              axis labels, respectively.
###
###      sides:  0:  Scale is -rng to +rng
###              1:  Scale is 0 to +rng
###             -1:  Scale is -rng to 0
###
###  3.  plt.map(dat, coord, labels, cols, legends, norm, sp, …):  Plots points
###      on a projection of Antarctica.
###
###      dat:  The data to be plotted.
###
###      coord:  Lon/Lat coordinates of the points.  Longitude should be in column
###              1 and latitude in column 2.
###
###      cols:  Vector of colors.  Default is NA, which results in automatic
###             generation of a color scale.
###
###      legends:  Default is NA, which results in automatic generation of a legend
###                for the color scale.  If manual legends are desired, legends
###                requires a list of:
###                [[1]]  Legend location (i.e., “bottomleft”)
###                [[2]]  Vector of legend text
###                [[3]]  Value for pch
###                [[4]]  Vector of colors corresponding to legend text
###                [[5]]  cex value for legend text and points
###                [[6]]  Background color
###
###      norm:  If TRUE, normalizes the plotted data to a scale of +/- 1.
###
###      sp:  Size parameter for plotted points.
###
###      … :  Parameters to be passed to set.colors and set.legend for automatic
###             color map and legend generation.
###

set.colors=function(rng, divs, hlite, hcol, sides) {

###  Make color scale
cold.colors=matrix(nrow=divs); hot.colors=matrix(nrow=divs)
for(i in 1:divs) {
lum.ratio=25*(i/divs)
c1=ifelse(lum.ratio<15, 255-(17*lum.ratio), 0)
c2=ifelse(lum.ratio<15, 255, 255-(10*(lum.ratio-15)))
h1=ifelse(255-(12*lum.ratio)>0, 255-(12*lum.ratio), 0)
h2=ifelse(255-(25*lum.ratio)>0, 255-(25*lum.ratio), 0)
cold.colors[divs+1-i]=rgb(red=c1, green=c1, blue=c2, maxColorValue=255)
hot.colors[i]=rgb(red=c2, green=h1, blue=h2, maxColorValue=255)
}

###  Define the zero point color.  If few enough divs for a space between each color,
###  then use gray for the zero point color.  Otherwise, use white.
neut.col=ifelse(divs<11, “#EEEEEE”, “#FFFFFF”)

###  Combine the colors into the color scale vector
if(sides==0) {
mapcolors=c(“darkmagenta”, as.vector(cold.colors), neut.col, as.vector(hot.colors), “black”) }
if(sides==-1) {
mapcolors=c(“darkmagenta”, as.vector(cold.colors), neut.col) }
if(sides==1) {
mapcolors=c(neut.col, as.vector(hot.colors), “black”) }

###  If hlite is used, replace appropriate portions of the color scale with hcol
if(!is.na(hlite[1])) {
hlite[2]=ifelse(is.na(hlite[2]), hlite[1], hlite[2])
start.pos=ifelse(sides==0, (divs+2), 1)
start.pos=ifelse(sides==-1, length(mapcolors), start.pos)
highlmin=start.pos+(round((divs+1)*(min(hlite[1], hlite[2])/(rng))))
highlmin=ifelse(highlmin<1, 1, highlmin)
highlmax=start.pos+(round((divs+1)*(max(hlite[1], hlite[2])/(rng))))
highlmax=ifelse(highlmax>length(mapcolors), length(mapcolors), highlmax)
mapcolors[highlmin:highlmax]=hcol
}

###  Return the color scale vector
mapcolors
}

set.legend=function(x, y, sl, st, spl, h, hcol, rng, ml, mbx, mxs, mc, divs, s.main, s.sub, s.axis, sides) {

###  Draw the box for the legend
rect(xleft=x, xright=x+(.6*sl), ytop=y, ybottom=y-(.080*sl), col=”white”)

###  If sides=0, plot the neutral color (exact center of the color scale)
if(sides==0) {points(x=x+(.3*sl), y=y-(.023*sl), col=mc[2+divs], pch=15, cex=1.3*spl)}

###  Plot the rest of the color scale
if(sides==-1 | sides==1) {
for(i in 1:divs) {
points(x=x+.1*sl+(.4/(divs+1))*i*sl, y=y-(.023*sl), col=mc[1+i], pch=15, cex=1.3*spl)
}
}
if(sides==0) {
for(i in 1:divs) {
points(x=x+(.3*sl)-(.2/(divs+1))*i*sl, y=y-(.023*sl), col=mc[divs+2-i], pch=15, cex=1.3*spl)
points(x=x+(.3*sl)+(.2/(divs+1))*i*sl, y=y-(.023*sl), col=mc[divs+2+i], pch=15, cex=1.3*spl)
}
}

###  Plot the -limit and +limit colors of the color scale
points(x=x+(.1*sl), y=y-(.023*sl), col=mc[1], pch=15, cex=1.3*spl)
points(x=x+(.5*sl), y=y-(.023*sl), col=mc[length(mc)], pch=15, cex=1.3*spl)

###  Plot the units, -limit label, and +limit label
if(sides==0) {lmt1=-rng;lmt2=rng}
if(sides==-1) {lmt1=-rng; lmt2=0}
if(sides==1) {lmt1=0;lmt2=rng}
text(x=x+(.3*sl), y=y-(.055*sl), ml[1], cex=st)
text(x=x+(.085*sl), y=y-(.055*sl), lmt1, cex=st)
text(x=x+(.51*sl), y=y-(.055*sl), lmt2, cex=st)

###  If hlite is used, define the text describing the highlighted range
if(!is.na(h[2])) {
ltext=paste(min(h[1],h[2]), “to”, max(h[1], h[2]))
ltext=ifelse(min(h[1],h[2])<=lmt1, paste(max(h[1],h[2]), “and Less”), ltext)
ltext=ifelse(max(h[1],h[2])>=lmt2, paste(min(h[1],h[2]), “and Greater”), ltext)
}
rtext=ifelse(!is.na(h[2]), ltext, h[1])
xtext=ifelse(!is.na(h[1]), paste(rtext, ml[1], “Highlighted”), NA)

###  Print the title, subtitle, and highlighted range text
title(main=ml[2], sub=ml[3], xlab=xtext, cex.main=s.main, cex.sub=s.sub, cex.axis=s.axis)

###  Draw box and/or axes?
if(mbx==TRUE) {box()}
if(mxs==TRUE) {map.axes()}
}

map.lbls=c(“Deg C/Decade”, “”, “”)

plt.map=function(dat, coord=sat.coord, rng=.5, divs=1000, labels=map.lbls, norm=FALSE, hlite=NA, hcol=”#99FF33″, mbx=TRUE, mxs=FALSE, x=-.3, y=-.35, sl=1, st=1, sp=1, spl=1, s.main=1, s.sub=1, s.axis=1, cols=NA, legends=NA, sides=0) {

library(mapproj)

###  Set up Antarctic projection and lat/lon lines
map(“world”, proj=”azequidistant”, orient=c(-90, 0, 0), ylim=c(-90, -60))
a=seq(-180, 150, 30)
for(i in 1:length(a)) {lines(mapproject(list(x=rep(a[i], 181), y=c(-90:90))), col=4, lty=2)}
a=c(-180,0)
for(i in 1:length(a)) {lines(mapproject(list(x=rep(a[i], 181), y=c(-90:90))), col=4, lty=1)}
a=seq(-50,-80,-10)
for(i in 1:4) {lines(mapproject(list(y=rep(a[i], 72), x=seq(-177.5, 177.5,5))), col=4, lty=2)}

###  If no user provided color scale, make the color scale
if(is.list(cols)) {mapcolors=cols} else {mapcolors=set.colors(rng, divs, hlite, hcol, sides)}

###  Normalize the data to a maximum of +/- 1?
nm=ifelse(norm==TRUE, max(abs(as.vector(dat))), 1)

###  Assign the color scale
if(!is.list(cols)) {
start.pos=ifelse(sides==0, (divs+2), 1)
start.pos=ifelse(sides==-1, length(mapcolors), start.pos)
temp=start.pos+(round((divs+1)*(dat/(rng*nm))))
temp=ifelse(temp<1, 1, temp)
temp=ifelse(temp>length(mapcolors), mapcolors[length(mapcolors)], mapcolors[temp])
} else {temp=cols}

###  Plot points
xy=mapproject(coord[, 1], coord[, 2])
points(xy$x, xy$y, pch=15, cex=sp, col=temp)

###  Overlay map
map(“world”, proj=”azequidistant”, orient=c(-90, 0, 0), ylim=c(-90, -60), fill=FALSE, add=TRUE)

###  Set up legend and labels
if(!is.list(legends)) {
legends=set.legend(x, y, sl, st, spl, h=hlite, hcol, rng, ml=labels, mbx, mxs, mc=mapcolors, divs, s.main, s.sub, s.axis, sides)
} else {
legend(legends[[1]], legends[[2]], pch=legends[[3]], col=legends[[4]], cex=legends[[5]], bg=legends[[6]])
title(main=labels[1], sub=labels[2])
}
}

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

calc.anom=function(dat) {

###  A present from Steve
for(i in 1:12) { sequ = seq(i, nrow(dat), 12)
dat[sequ, ] = scale(dat[sequ, ], scale=F) }
dat
}

ssq=function(x) {sum(x^2)}

regem.IPCA=function(X, maxiter=5000, s.tol=.05, f.tol=.01, regpar=3, method=”default”, p.info=”Predictor Matrix”) {

###  Define regem, which contains the output of each iteration
regem=list()

###  Regpar cannot be greater than full rank
regpar=ifelse(regpar>ncol(X), ncol(X), regpar)

###  Find matrix dimensions and missing values
n=nrow(X);p=ncol(X)
indmis=is.na(X)
nmis=sum(indmis)
dofC=n-1
dofS=dofC-min(regpar, n-1, p)

###  Center and infill NAs with zeros (prep for initial infill)
X=scale(X, scale=FALSE)
M=attributes(X)$”scaled:center”
X[indmis]=0

###  Do initial infill
init.regem=regem.sub(X, 1, 1, s.tol, indmis, nmis, p, dofC, dofS, M, “terse”)

###  Extract infilled matrix
X=init.regem$X

###  Remove old scaling/centering values
attr(X, “scaled:center”)=NULL
attr(X, “scaled:scale”)=NULL

###  Recenter
X=scale(X, scale=FALSE)
M=attributes(X)$”scaled:center”

###  Define loop start
for(regp in 1:regpar) {

###  If on final iteration, use the tighter tolerance.
tol=ifelse(regp<regpar, s.tol, f.tol)

###  Impute using rank “regp” (slow approach to full rank)
regem[[regp]]=regem.sub(X, regp, maxiter, tol, indmis, nmis, p, dofC, dofS, M, method)

###  Extract results
subreg=regem[[regp]]

###  Extract final iteration (method “terse” only returns the final iteration)
if(method==”default”) {
L=length(subreg)
subreg=subreg[[L]]
}

###  Extract imputed matrix and remove old scaling/centering values
X=subreg$X
attr(X, “scaled:center”)=NULL
attr(X, “scaled:scale”)=NULL

###  Recenter to prepare for next loop
X=scale(X, scale=FALSE)
M=attributes(X)$”scaled:center”

###  Print status
print(p.info, quote=FALSE)
print(“———————–“, quote=FALSE)
print(paste(“Imputation at regpar=”, regp, ” complete; “, subreg$info, sep=””), quote=FALSE)
print(paste(“dXmis=”, subreg$dXmis, sep=””), quote=FALSE); print(“”, quote=FALSE)
}
###  Return all iterations
regem
}

regem.sub=function(X, regp, maxiter, tol, indmis, nmis, p, dofC, dofS, M, method) {

###  Define subreg, which contains the output of each iteration
subreg=list()

###  Set up loop
iter=0
dXmis=10^6

while(iter<=maxiter &  dXmis>tol) {

###  Increment
iter=iter+1

###  Define cov(X)
C=t(X) %*% X/dofC
CovRes=array(0, dim=c(p, p))

###  Define sd
D=sqrt(diag(C))

###  Normalize to sd=1
X=X %*% diag(1/D)
C=diag(1/D )%*% C %*% diag(1/D)

###  Find the singular value decomposition of the input matrix
###  using regp to truncate the number of right/left singular
###  vectors to be computed and calculate the truncated rank matrix.
svd.X= svd(X, nu=regp, nv=regp)
if(regp>1) {
Xmis= svd.X$u %*% diag(svd.X$d[1:regp]) %*% t(svd.X$v)
}
if(regp==1) {
Xmis= as.matrix(svd.X$u) %*% as.matrix(svd.X$d[1]) %*% t(svd.X$v)
}

###  Rescale to original scaling
X=X1a=scale(X, center=FALSE, scale=1/D)
Xmis=scale(Xmis, center=FALSE, scale=1/D)
C= diag(D) %*% C %*% diag(D)
CovRes=diag(D) %*% CovRes %*% diag(D)

###  Find iteration imputation error
dXmis=sqrt(ssq(Xmis[indmis]-X[indmis]))/sqrt(nmis)

###  Update X with imputed values
X[indmis]=Xmis[indmis]

###  Recenter
X=scale(X, scale=FALSE)
Mup=attributes(X)$”scaled:center”
M=M+Mup
C=t(X) %*% X/dofC

###  Store the iteration results and end the iteration
subreg[[iter]]=list(X=X, dXmis=dXmis)
}

if(dXmis>tol) {
conv.info=paste(iter, “iterations without convergence.”) } else {
conv.info=paste(iter, “iterations to convergence.”)
}

###  Store the SVD of X after infilling and store the solution matrix
N=length(subreg)
subreg[[N]]$svd=svd.X
subreg[[N]]$Xmis=Xmis
subreg[[N]]$info=conv.info

###  Method “terse” saves only the final iteration for each value of regp
if(method==”terse”) {
regem=subreg[[N]] } else {
regem=subreg
}
###  Return the full set of iteration information
regem
}

get.PCs=function(dat, PC.max=3, st=1982) {

###  Calculate the singular value decomposition of the input matrix
dat=t(dat)
dat.svd=svd(dat)

###  Multiply the right singular vectors through and store as a time series
PCs=ts(t((dat.svd$d)*t(dat.svd$v)), start=st, freq=12)

###  Print explained variance
var.tot=0
print(“PC Selection”, quote=FALSE); print(“============”, quote=FALSE)
print(“Explained variance:”, quote=FALSE); print(“”, quote=FALSE)
for(i in 1:PC.max) {
print(paste(“PC “, i, “:  “, round((dat.svd$d[i]^2/ssq(dat.svd$d))*100,2), “%”, sep=””), quote=FALSE)
var.tot=var.tot+(dat.svd$d[i]^2/ssq(dat.svd$d))*100
}
print(“”, quote=FALSE); print(paste(“Total explained variance:  “, round(var.tot,2), “%”, sep=””), quote=FALSE); print(“”, quote=FALSE)

###  Return the desired number of PCs
PCs=list(PCs[, 1:PC.max], dat.svd$u[, 1:PC.max], dat.svd$v, dat.svd$d)
PCs
}

get.recon=function(dat.predict, dat.cov, PC.max) {

###  Multiply the left singular vectors (spatial weights) through
if(PC.max==1) {

###  If only 1 PC, dat.cov is a vector
recon=array(0, dim=c(600, length(dat.cov)))
for(i in 1:length(dat.cov)){

###  Convert to time series and store
recon[, i]=ts(dat.predict*dat.cov[i], start=1957, freq=12)
} } else {

###  If more than 1 PC, dat.cov is a matrix
recon=array(0, dim=c(600, nrow(dat.cov)))
for(i in 1:nrow(dat.cov)) {

###  Convert to time series
temp.rec=ts(t(t(dat.predict)*(dat.cov[i, ])), start=1957, freq=12)

###  Sum the individual PCs at each location and store
recon[,i]=ts(rowSums(temp.rec), start=1957, freq=12)
}
}
###  Return the reconstruction
recon
}

stat.trends=function(dat, st, en) {
dat=window(dat, st, en)

###  Store trends as deg C per decade
temp=lm(dat~I(time(dat)))[[1]][2, ]*10
temp
}

info=c(“”, “Elements in returned list:”,”=========================”, “[[1]]:  Reconstruction time series”,”[[2]]:  Reconstruction overall monthly means time series”,”[[3]]:  Imputed PCs”,”[[4]]:  Final RegEM iteration for each truncation setting of the predictor matrix”,”[[5]][[1]]:  Ground station and PCs matrix (input to RegEM)”,”[[5]][[2]]:  Left singular vectors (1 – PC.max)”,”[[5]][[3]]:  Right singular vectors (all)”,”[[5]][[4]]:  Eigenvalues (all)”, “[[6]]:  Linear trends (1957-2006)”, “[[7]]:  Linear trends (1965-2006)”, “[[8]]:  Linear trends (1982-2006)”, “[[9]]:  Linear trends (1982-1994.5)”, “[[10]]:  Linear trends (1994.5-2006)”)

do.recon=function(dat=avhrr.anom, base=reader.steig.anom, PC.max=3, maxiter=4000, s.tol=.05, f.tol=.01, regpar.dat=3, regpar.base=3, method=”terse”, base.st=1957, dat.st=1982, rec.en=c(2006,12), infor=info, tr=TRUE, split=FALSE) {

###  Get the PCs and SVD
PCs=get.PCs(dat, PC.max, dat.st)

###  Extract the PCs for imputation
PC.dat=PCs[[1]]

###  Extract the spatial weights (left eigenvector)
PC.cov=PCs[[2]]

###  Impute the predictor matrix
regpar1=ifelse(regpar.base>ncol(base), ncol(base), regpar.base)
gnd.regem=regem.IPCA(base, maxiter, s.tol, f.tol, regpar1, method=”terse”)

###  Extract the result, remove centering/scaling notes, store as a time series
X=gnd.regem[[regpar1]]$X
attr(X, “scaled:center”)=NULL
attr(X, “scaled:scale”)=NULL
X=ts(X, start=base.st, freq=12)

###  Set up storage for the imputed PCs
imputed.pcs=matrix(nrow=nrow(base), ncol=PC.max)

###  If split=TRUE, impute each PC independently
if(split==TRUE) {
for(i in 1:PC.max) {

###  Place PC number “i” alongside the predictor matrix
Xi=ts.union(X, PC.dat[, i])

###  Impute the PC and store
regpar1=ifelse(regpar.dat>ncol(Xi), ncol(Xi), regpar.dat)
pc.regem=regem.IPCA(Xi, maxiter, s.tol, f.tol, regpar1, method=”terse”, p.info=paste(“PC Imputation; PC”, i))
imputed.pcs[,i]=pc.regem[[PC.max]]$X[, ncol(base)+1]
}
}

###  If split=FALSE, impute all PCs together
if(split!=TRUE) {

###  Place all PCs alongside the predictor matrix
Xi=ts.union(X, PC.dat)

###  Impute the PCs and store
regpar1=ifelse(regpar.dat>ncol(Xi), ncol(Xi), regpar.dat)
pc.regem=regem.IPCA(Xi, maxiter, s.tol, f.tol, regpar1, method=”terse”, p.info=”PC Imputation”)
imputed.pcs=pc.regem[[regpar1]]$X[, -(1:ncol(X))]
}

###  Multiply the left eigenvector through to perform the recon
recon=ts(get.recon(imputed.pcs, PC.cov, PC.max), start=base.st, freq=12)

###  Calculate standard trends
trend.all=0; trend.65=0; trend.sat=0; trend.82=0; trend.94=0
if(tr==TRUE) {
trend.all=stat.trends(recon, 1957, c(2006,12))
trend.65=stat.trends(recon, 1965, c(2006,12))
trend.sat=stat.trends(recon, 1982, c(2006,12))
trend.82=stat.trends(recon, 1982, c(1994,6))
trend.94=stat.trends(recon, c(1994,7), c(2006,12))
}

###  Calculate continent-wide average
recon.means=ts(rowMeans(recon), start=base.st, freq=12)

###  Collate the recon, trends, SVD information, and RegEM information
all.recon=list(recon, recon.means, imputed.pcs, gnd.regem, PCs, trend.all, trend.65, trend.sat, trend.82, trend.94)

###  Print the list description and return the list
print.table(info)
all.recon
}

###  Load files

lon=read.table(“tir_lons.txt”)
lat=read.table(“tir_lats.txt”)
sat.coord = cbind(lon,lat)
load(“reader.data.04-04-2009.RData”)
load(“reader.steig.04-12-2009.RData”)
load(“reader.jeff.RData”)
grid=scan(“temp.dat”, n= -1)
avhrr=ts(t(array(grid, dim=c(5509, 300))), start=1982, freq=12); rm(grid)
grid=scan(“ant_recon.txt”, n= -1)
ant.recon=ts(t(array(grid, dim=c(5509, 600))), start=1957, freq=12); rm(grid)

###  Calculate anomalies
reader.anom=calc.anom(reader.data); reader.steig.anom=calc.anom(reader.steig); avhrr.anom=calc.anom(avhrr)
rm(avhrr, reader.data, reader.steig, lon, lat)

###  Get AVHRR and main TIR recon trends
#  Use to calculate:

# m.trends=list()
# a.trends=list()

# m.trends[[6]]=stat.trends(ant.recon, 1957, c(2006,12))
# m.trends[[7]]=stat.trends(ant.recon, 1965, c(2006,12))
# m.trends[[8]]=stat.trends(ant.recon, 1982, c(2006,12))
# m.trends[[9]]=stat.trends(ant.recon, 1982, c(1994,6))
# m.trends[[10]]=stat.trends(ant.recon, c(1994,7), c(2006,12))

# a.trends[[6]]=0
# a.trends[[7]]=0
# a.trends[[8]]=stat.trends(avhrr.anom, 1982, c(2006,12))
# a.trends[[9]]=stat.trends(avhrr.anom, 1982, c(1994,6))
# a.trends[[10]]=stat.trends(avhrr.anom, c(1994,7), c(2006,12))

#  Alternately load saved trends:

#load(“avhrr.trends.RData”)
#load(“TIR.trends.RData”)

#m.means=ts(rowMeans(ant.recon), start=1957, freq=12)
#a.means=ts(rowMeans(avhrr.anom), start=1982, freq=12)

###  Do a 3-PC, 3-regpar recon:

#recon=do.recon()