## This script contains the analysis code used precisely by than manuscript

library(fields)
library(maps)
library(rgl)
library(abind)
library(ncdf)

cdd1=read.csv(file="CDDNA.csv", header=T)[,-1]
hdd1=read.csv(file="HDDNA.csv", header=T)[,-1]
coord=read.csv(file="coordNA.csv", header=T)[,-1]
cities=read.csv(file="cities.csv", header=T)

nc=open.ncdf("orog_CRCM.nc")
orog=get.var.ncdf(nc, "orog")
oro=orog[,1:85]

lon=matrix(0,140,115)
for (i in 1:115) {
  for (j in 1:140) {
    count=(j-1)*115+i
    lon[j,i]=coord[count,1]
  }
}

lat=matrix(0,140,115)
for (i in 1:115) {
  for (j in 1:140) {
    count=(j-1)*115+i
    lat[j,i]=coord[count,2]
  }
}

loc=cbind(c(lon[,1:85]),c(lat[,1:85]))

rowkeep=seq(1,140,by=2)
colkeep=seq(1,85,by=2)
locu=cbind(c(lon[rowkeep,colkeep]),c(lat[rowkeep,colkeep]))

cdd=array(0,dim=c(140,115,31))
for (y in 1:31) {
  for (i in 1:115) {
    for (j in 1:140) {
      print(y + i/1000)
      count=(i-1)*140+j
      cdd[j,i,y]=cdd1[count,y]
    }
  }
}

hdd=array(0,dim=c(140,115,30))
for (y in 1:30) {
  for (i in 1:115) {
    for (j in 1:140) {
      print(y + i/1000)
      count=(i-1)*140+j
      hdd[j,i,y]=hdd1[count,y]
    }
  }
}


######################################################
## Complete Spatial Analysis of Trends in CDD
## 1. Compute initial estimate of Trend, then detrend data to year 2038

seq1=seq(1,31,by=1)
cdd.mean=apply(cdd[1:140,1:85,],c(1,2),mean)

cdd.alpha.twiddle=matrix(0,140,85)
cdd.beta.twiddle=matrix(0,140,85)
cdd.beta.se=matrix(0,140,85)
cdd.sd=matrix(0,140,85)
for (i in 1:140) {
  for (j in 1:85) {
    print(i+j/1000)
    cdd.alpha.twiddle[i,j]=summary(lm(cdd[i,j,]~seq1))$coeff[1,1]
    cdd.beta.twiddle[i,j]=summary(lm(cdd[i,j,]~seq1))$coeff[2,1]
    cdd.beta.se[i,j]=summary(lm(cdd[i,j,]~seq1))$coeff[2,2]
    cdd.sd[i,j]=summary(lm(cdd[i,j,]~seq1))$sigma
  }
}

## Empirical Estimates of Mean, SD, and Trend.  This is called "Emp1.eps"
par(mfrow=c(2,2))
quilt.plot(c(loc[,1]),c(loc[,2]),c(cdd.mean),nx=100,ny=73, main="Empirical Estimate of Mean",xlab="Lon",ylab="Lat",cex.lab=1.3,cex.main=1.3,cex.axis=1.3,cex.index=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

quilt.plot(c(loc[,1]),c(loc[,2]),c(cdd.beta.twiddle),nx=100,ny=73, main="Empirical Estimate of Trend",xlab="Lon",ylab="",cex.lab=1.3,cex.main=1.3,cex.axis=1.3,cex.index=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

quilt.plot(c(loc[,1]),c(loc[,2]),c(cdd.sd),nx=100,ny=73,, main="Empirical Estimate of Std. Dev.",xlab="Lon",ylab="Lat",cex.lab=1.3,cex.main=1.3,cex.axis=1.3,cex.index=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
dev.off()

## Which grid cells can we even estimate trend surfaces for?  Count number of zeroes in data.
zero=matrix(0,140,85)
for (i in 1:140) {
  for (j in 1:85) {
    zero[i,j]=sum(cdd[i,j,]>0)
  }
}

## Rescale data.  This will give the value "NA" if there are fewer than 25 non-zero pojected CDD totals.  
cdd.rescale=cdd[,1:85,]-cdd[,1:85,]
for (i in 1:140) {
  for (j in 1:85) {
    for (k in 1:31) {
      print(i+j/1000)
      if (zero[i,j]>24) {cdd.rescale[i,j,k]=(cdd[i,j,k] - (k-16)*cdd.beta.twiddle[i,j] - cdd.mean[i,j])/cdd.sd[i,j]} else {cdd.rescale[i,j,k] <- NA}
    }
  }
}

## A sample residual.  This is called "Emp2.png".  It's important to save as PNG to flatten the image.
par(mfrow=c(1,2))
k=3
quilt.plot(c(lon[,1:85]),c(lat[,1:85]),c(cdd.rescale[,,k]), main="Residuals from Standardized Data in Year 2040",zlim=c(-4,4),xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

qqnorm(cdd.rescale[,,],main="QQ Plot of All Std. Residuals")
qqline(cdd.rescale[,,])
dev.off()

## Remove some unneeded objects to free up memory
rm(cdd1)
rm(hdd1)
rm(zero)



## Without orography
cmu2=Tps(locu,c(cdd.mean[rowkeep,colkeep]), lon.lat=TRUE)
csdu2=Tps(locu,c(cdd.sd[rowkeep,colkeep]), lon.lat=TRUE)
predcmu2 <- predict(cmu2, loc)
predcsd2 <- predict(csdu2, loc)

rm(cmu2)
rm(csdu2)

## With orography
cmu=Tps(locu,c(cdd.mean[rowkeep,colkeep]),Z=c(oro[rowkeep,colkeep]), lon.lat=TRUE)
csdu=Tps(locu,c(cdd.sd[rowkeep,colkeep]),Z=c(oro[rowkeep,colkeep]), lon.lat=TRUE)
predcmu <- predict(cmu, loc, Z=c(oro))
predcsd <- predict(csdu, loc, Z=c(oro))

## Look at fit in mean
par(mfrow=c(1,2))
quilt.plot(loc,c(cdd.mean),main="Raw Mean CDD data",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,6500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predcmu,main="Thin Plate Spline fit to Mean with Orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,6500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(cdd.mean),main="Raw Mean CDD data",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,6500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predcmu2,main="Thin Plate Spline fit to Mean without Orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,6500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(cdd.mean)-predcmu2,main="Residuals in Mean fit without orography",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(cdd.mean)-predcmu,main="Residualsin Mean fit with orography",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

## This is a Publication Image
par(mfrow=c(1,2))
quilt.plot(loc,predcmu,main="Thin Plate Spline fit to Mean",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,6500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(cdd.mean)-predcmu,main="Residuals",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

persp3d(x=seq(1,140,by=1),y=seq(1,85,by=1),z=matrix(c(cdd.mean)-predcmu,nrow=140,ncol=85),
zlab="Residual",xlab="Lon",ylab="Lat",col="grey")

## Look at fit in SD
par(mfrow=c(1,2))
quilt.plot(loc,c(cdd.sd),main="Raw SD CDD data",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predcsd2,main="Thin Plate Spline fit to SD without orography",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(cdd.sd),main="Raw SD CDD data",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predcsd,main="Thin Plate Spline fit to SD with orography",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(cdd.sd)-predcsd2,main="Residuals in SD fit without orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-150,150))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(cdd.sd)-predcsd,main="Residuals in SD fit with orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-150,150))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

## Reduction in MSE in Mean TPS fit to CDD
mean((c(cdd.mean)-predcmu2)^2)
mean((c(cdd.mean)-predcmu)^2)
(mean((c(cdd.mean)-predcmu2)^2)-mean((c(cdd.mean)-predcmu)^2))/(mean((c(cdd.mean)-predcmu2)^2))

## This is a Publication Image
par(mfrow=c(1,2))
quilt.plot(loc,predcsd,main="Thin Plate Spline fit to SD",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(cdd.sd)-predcsd,main="Residuals",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-100,100))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

persp3d(x=seq(1,140,by=1),y=seq(1,85,by=1),z=matrix(c(cdd.sd)-predcsd,nrow=140,ncol=85),
zlab="Residual",xlab="Lon",ylab="Lat",col="grey")

## Single Publication Image, "TPSCDDFits.eps"
par(mfrow=c(2,2))
quilt.plot(loc,predcmu,main="Thin Plate Spline fit to Mean of CDD",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,6500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(cdd.mean)-predcmu,main="Residuals",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

quilt.plot(loc,predcsd,main="Thin Plate Spline fit to SD of CDD",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(cdd.sd)-predcsd,main="Residuals",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-100,100))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)


#######################################
## The grand fit!!!

zu=(c(cdd[rowkeep,colkeep,1]) - cmu$fitted.values) / csdu$fitted.values
fit = Krig(locu, zu, Z=c(oro[rowkeep,colkeep])) 
pred3 <- predict(fit, loc, drop.Z=TRUE) 
surface = pred3*predcsd + predcmu

cout=matrix(as.vector(surface),nrow=140,ncol=85)

for (i in 2:31) {
  print(i)
  zu=(c(cdd[rowkeep,colkeep,i]) - cmu$fitted.values) / csdu$fitted.values
  fit <- Krig(locu, zu, Z=c(oro[rowkeep,colkeep])) 
  pred3 <- predict(fit, loc, drop.Z=TRUE) 
  surface = pred3*predcsd + predcmu
  p4=matrix(as.vector(surface),nrow=140,ncol=85)
  cout=abind(cout,p4,along=3)
}

slope.cdd=matrix(0,140,85)
se.slope.cdd=matrix(0,140,85)
seq1=seq(1,31,by=1)

for (i in 1:140) {
  for (j in 1:85) {
    print(i+j/1000)
    slope.cdd[i,j]=summary(lm(cout[i,j,]~seq1))$coeff[2,1]
    se.slope.cdd[i,j]=summary(lm(cout[i,j,]~seq1))$coeff[2,2]
  }
}

quilt.plot(loc,surface)
par(mfrow=c(1,2))
z=c(cdd[1:140,1:85,31])
quilt.plot(loc,z-cout[,,20],zlim=c(-max(abs(z-cout[,,20])),max(abs(z-cout[,,20]))))
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
hist(z-surface, main="Residuals",xlab="",ylab="")
persp3d(x=seq(1,140,by=1),y=seq(1,85,by=1),z=matrix(z-surface,nrow=140,ncol=85),
zlab="Residual",xlab="Lat",ylab="Lon",col="grey")
mean((z-surface)^2)

par(mfrow=c(2,2)) ## This is "CDDTrend.eps"
quilt.plot(loc,c(slope.cdd),nx=100,ny=85,main="Estimated Trends in CDD",xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(se.slope.cdd),nx=100,ny=85,main="Estimated Std. Err. of Trends in CDD",xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

quilt.plot(loc,c(slope.cdd)-2*c(se.slope.cdd),nx=100,ny=85,zlim=c(-15,45),main="Lower 95% Confidence Surface",xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(slope.cdd)+2*c(se.slope.cdd),nx=100,ny=85,zlim=c(-15,45),main="Upper 95% Confidence Surface",xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)


######################################################
## Complete Spatial Analysis of Trends in HDD

seq1=seq(1,30,by=1)
hdd.mean=apply(hdd[1:140,1:85,],c(1,2),mean)

hdd.alpha.twiddle=matrix(0,140,85)
hdd.beta.twiddle=matrix(0,140,85)
hdd.beta.se=matrix(0,140,85)
hdd.sd=matrix(0,140,85)
for (i in 1:140) {
  for (j in 1:85) {
    print(i+j/1000)
    hdd.alpha.twiddle[i,j]=summary(lm(hdd[i,j,]~seq1))$coeff[1,1]
    hdd.beta.twiddle[i,j]=summary(lm(hdd[i,j,]~seq1))$coeff[2,1]
    hdd.beta.se[i,j]=summary(lm(hdd[i,j,]~seq1))$coeff[2,2]
    hdd.sd[i,j]=summary(lm(hdd[i,j,]~seq1))$sigma
  }
}

## Empirical Estimates of Mean, SD, and Trend

par(mfrow=c(1,2))
quilt.plot(c(loc[,1]),c(loc[,2]),c(hdd.mean),nx=100,ny=73, main="Empirical Estimate of Mean",xlab="Lon",ylab="Lat",cex.lab=1.3,cex.main=1.3,cex.axis=1.3,cex.index=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(c(loc[,1]),c(loc[,2]),c(hdd.sd),nx=100,ny=73,, main="Empirical Estimate of Std. Dev.",xlab="Lon",ylab="Lat",cex.lab=1.3,cex.main=1.3,cex.axis=1.3,cex.index=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
dev.off()

quilt.plot(c(loc[,1]),c(loc[,2]),c(hdd.beta.twiddle),nx=100,ny=73, main="Empirical Estimate of Trend",xlab="Lon",ylab="Lat",cex.lab=1.3,cex.main=1.3,cex.axis=1.3,cex.index=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

## Which grid cells can we even estimate trend surfaces for?  Count number of zeroes in data.
hzero=matrix(0,140,85)
for (i in 1:140) {
  for (j in 1:85) {
    hzero[i,j]=sum(hdd[i,j,]>0)
  }
}

## Rescale data.  This will give the value "NA" if there are fewer than 25 non-zero pojected CDD totals.  
hdd.rescale=hdd[,1:85,]-hdd[,1:85,]
for (i in 1:140) {
  for (j in 1:85) {
    for (k in 1:30) {
      print(i+j/1000)
      if (hzero[i,j]>24) {hdd.rescale[i,j,k]=(hdd[i,j,k] - (k-16)*hdd.beta.twiddle[i,j] - hdd.mean[i,j])/hdd.sd[i,j]} else {hdd.rescale[i,j,k] <- NA}
    }
  }
}

par(mfrow=c(1,2))
k=3
quilt.plot(c(lon[,1:85]),c(lat[,1:85]),c(hdd.rescale[,,k]), main="Residuals from Standardized Data in Year 2040",zlim=c(-4,4),xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

qqnorm(hdd.rescale[,,],main="QQ Plot of All Std. Residuals")
qqline(hdd.rescale[,,])
dev.off()

## Without orography
hmu2=Tps(locu,c(hdd.mean[rowkeep,colkeep]), lon.lat=TRUE)
hsdu2=Tps(locu,c(hdd.sd[rowkeep,colkeep]), lon.lat=TRUE)
predhmu2 <- predict(hmu2, loc)
predhsd2 <- predict(hsdu2, loc)

## With orography
hmu=Tps(locu,c(hdd.mean[rowkeep,colkeep]),Z=c(oro[rowkeep,colkeep]), lon.lat=TRUE)
hsdu=Tps(locu,c(hdd.sd[rowkeep,colkeep]),Z=c(oro[rowkeep,colkeep]), lon.lat=TRUE)
predhmu <- predict(hmu, loc, Z=c(oro))
predhsd <- predict(hsdu, loc, Z=c(oro))

## Look at fit in mean
par(mfrow=c(1,2))
quilt.plot(loc,c(hdd.mean),main="Raw Mean HDD data",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,9500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predhmu,main="Thin Plate Spline fit to Mean with Orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,9500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(hdd.mean),main="Raw Mean HDD data",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,9500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predhmu2,main="Thin Plate Spline fit to Mean without Orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,9500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(hdd.mean)-predhmu2,main="Residuals in Mean fit without orography",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(hdd.mean)-predhmu,main="Residuals in Mean fit with orography",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

## Reduction in MSE in Mean TPS fit to HDD
mean((c(hdd.mean)-predhmu2)^2)
mean((c(hdd.mean)-predhmu)^2)
(mean((c(hdd.mean)-predhmu2)^2)-mean((c(hdd.mean)-predhmu)^2))/(mean((c(hdd.mean)-predhmu2)^2))

## This is a Publication Image
par(mfrow=c(1,2))
quilt.plot(loc,predhmu,main="Thin Plate Spline fit to Mean with Orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,9500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(hdd.mean)-predhmu,main="Residuals in Mean fit with Orography",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

persp3d(x=seq(1,140,by=1),y=seq(1,85,by=1),z=matrix(c(hdd.mean)-predhmu,nrow=140,ncol=85),
zlab="Residual",xlab="Lon",ylab="Lat",col="grey")

## Look at fit in SD
par(mfrow=c(1,2))
quilt.plot(loc,c(hdd.sd),main="Raw SD HDD data",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predhsd2,main="Thin Plate Spline fit to SD without orography",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(hdd.sd),main="Raw SD HDD data",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,predhsd,main="Thin Plate Spline fit to SD with orography",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

par(mfrow=c(1,2))
quilt.plot(loc,c(hdd.sd)-predhsd2,main="Residuals in SD fit without orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-150,150))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(hdd.sd)-predhsd,main="Residuals in SD fit with orography",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-150,150))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

## This is a Publication Image
par(mfrow=c(1,2))
quilt.plot(loc,predhsd,main="Thin Plate Spline fit to SD with orography",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(hdd.sd)-predhsd,main="Residuals",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-100,100))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

persp3d(x=seq(1,140,by=1),y=seq(1,85,by=1),z=matrix(c(hdd.sd)-predhsd,nrow=140,ncol=85),
zlab="Residual",xlab="Lon",ylab="Lat",col="grey")


## This is "TPSHDdFits.eps"
par(mfrow=c(2,2))
quilt.plot(loc,predhsd,main="Thin Plate Spline fit to Mean of HDD",xlab="Lon",ylab="Lat",nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(hdd.sd)-predhsd,main="Residuals",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-100,100))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)

quilt.plot(loc,predhmu,main="Thin Plate Spline fit to SD of HDD",xlab="Lon",ylab="Lat",nx=100,ny=85,zlim=c(-500,9500))
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(hdd.mean)-predhmu,main="Residuals",xlab="Lon",ylab="Lat",zlim=c(-600,600),nx=100,ny=85)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-40), ylim=c(22, 60), lwd=2, add=T)


#######################################
## The grand fit!!!

zu=(c(hdd[rowkeep,colkeep,1]) - hmu$fitted.values) / hsdu$fitted.values
fit = Krig(locu, zu, Z=c(oro[rowkeep,colkeep])) 
pred3 <- predict(fit, loc, drop.Z=TRUE) 
surface1 = pred3*predhsd + predhmu

hout=matrix(as.vector(surface1),nrow=140,ncol=85)

for (i in 2:30) {
  print(i)
  zu=(c(hdd[rowkeep,colkeep,i]) - hmu$fitted.values) / hsdu$fitted.values
  fit <- Krig(locu, zu, Z=c(oro[rowkeep,colkeep])) 
  pred3 <- predict(fit, loc, drop.Z=TRUE) 
  surface = pred3*predhsd + predhmu
  p4=matrix(as.vector(surface),nrow=140,ncol=85)
  hout=abind(hout,p4,along=3)
}

image.plot(matrix(as.vector(surface),nrow=140,ncol=85)-matrix(as.vector(surface1),nrow=140,ncol=85),main="Change in HDD from 2038 to 2068")

slope.hdd=matrix(0,140,85)
se.slope.hdd=matrix(0,140,85)
seq1=seq(1,30,by=1)

for (i in 1:140) {
  for (j in 1:85) {
    print(i+j/1000)
    slope.hdd[i,j]=summary(lm(hout[i,j,]~seq1))$coeff[2,1]
    se.slope.hdd[i,j]=summary(lm(hout[i,j,]~seq1))$coeff[2,2]
  }
}

plot(seq1,hout[100,60,], type="l")

image.plot(slope.hdd)
par(mfrow=c(1,2))
image.plot(se.slope.cdd)
image.plot(se.slope.hdd)

quilt.plot(loc,surface)
par(mfrow=c(1,2))
z=c(hdd[1:140,1:85,30])
quilt.plot(loc,z-surface,zlim=c(-max(abs(z-surface)),max(abs(z-surface))))
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
hist(z-surface, main="Residuals",xlab="",ylab="")
persp3d(x=seq(1,140,by=1),y=seq(1,85,by=1),z=matrix(z-surface,nrow=140,ncol=85),
zlab="Residual",xlab="Lat",ylab="Lon",col="grey")
mean((z-pred3)^2)

par(mfrow=c(2,2)) ## This is "HDDTrend.eps"
quilt.plot(loc,c(slope.hdd),nx=100,ny=85,main="Estimated Trends in HDD")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(se.slope.hdd),nx=100,ny=85,main="Estimated Std. Err. of Trends in HDD")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)

quilt.plot(loc,c(slope.hdd)-2*c(se.slope.hdd),nx=100,ny=85,zlim=c(-40,15),main="Lower 95% Confidence Surface",xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
quilt.plot(loc,c(slope.hdd)+2*c(se.slope.hdd),nx=100,ny=85,zlim=c(-40,15),main="Upper 95% Confidence Surface",xlab="Lon",ylab="Lat")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)


##########################################
## This bit extracts the information for 30 North American cities
cities=read.csv(file="cities.csv", header=T)

quilt.plot(loc,c(slope.hdd),nx=100,ny=85,main="Estimated Trends")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
points(cities[,2],cities[,3],pch=19)

extract=matrix(0,30,2)
for (t in 1:30) {
  print(t)
  di=matrix(0,140,85)
  for (i in 1:140) {
    for (j in 1:85) {
      di[i,j]=sqrt((lon[i,j]-cities[t,2])^2+(lat[i,j]-cities[t,3])^2)
    }
  }
  extract[t,1]=which(di==min(di),arr.ind=TRUE)[1]
  extract[t,2]=which(di==min(di),arr.ind=TRUE)[2]
}

conf=1.98

cy=matrix(0,30,8)
colnames(cy)=c("CDD Trend","S.E.","Lower 95% Bound","Upper 95% Bound","HDD Trend","S.E.","Lower 95% Bound","Upper 95% Bound")
for (t in 1:30) {
  cy[t,1]=round(slope.cdd[extract[t,1],extract[t,2]],2)
  cy[t,2]=round(se.slope.cdd[extract[t,1],extract[t,2]],2)
  cy[t,3]=round(slope.cdd[extract[t,1],extract[t,2]]-conf*se.slope.cdd[extract[t,1],extract[t,2]],2)
  cy[t,4]=round(slope.cdd[extract[t,1],extract[t,2]]+conf*se.slope.cdd[extract[t,1],extract[t,2]],2)
  cy[t,5]=round(slope.hdd[extract[t,1],extract[t,2]],2)
  cy[t,6]=round(se.slope.hdd[extract[t,1],extract[t,2]],2)
  cy[t,7]=round(slope.hdd[extract[t,1],extract[t,2]]-conf*se.slope.hdd[extract[t,1],extract[t,2]],2)
  cy[t,8]=round(slope.hdd[extract[t,1],extract[t,2]]+conf*se.slope.hdd[extract[t,1],extract[t,2]],2)
}
data.frame(cities[,1],cy)

## The spatially smooth model gives average SE
mean(cy[,2])
mean(cy[,6])

## The pointwise, single grid cell approach gives average SE
mean(cdd.beta.se[extract[,1],extract[,2]])
mean(hdd.beta.se[extract[,1],extract[,2]])

(mean(cy[,6]) - mean(hdd.beta.se[extract[,1],extract[,2]]))/mean(hdd.beta.se[extract[,1],extract[,2]])

## All grid cells
mean(se.slope.cdd)
mean(cdd.beta.se[,])

mean(se.slope.hdd)
mean(hdd.beta.se[,])


wh=17
par(mfrow=c(1,3))
#quilt.plot(loc,c(slope.cdd),nx=100,ny=85,main="Locations of 30 Cities Traded on CME",col="white",add.legend="FALSE",xlab="Lon",ylab="Lat"
#,cex.lab=1.3,cex.main=1.3,cex.axis=1.3)
map(database="world", regions=c("usa","canada"), xlim=c(-140,-50), ylim=c(25, 60), lwd=1,cex.lab=1.5)
map.axes()
text(-138,63,lab="Locations of the 30 Cities with Traded Tremperature Derivatives")
points(cities[,2],cities[,3],pch=3,cex=1,lwd=3)

dev.off()

par(mfrow=c(1,2)) ## This is "LocsNYC.eps"
plot(cout[extract[wh,1],extract[wh,2],]~seq(2038,2068,by=1),ylab="CDD",xlab="Year",type="l",lwd=1.5,main="Trend in New York CDDs",
cex.lab=1.5,cex.main=1.5,cex.axis=1.5,cex=1.3)
points(cout[extract[wh,1],extract[wh,2],]~seq(2038,2068,by=1),pch=19)
points(cdd[extract[wh,1],extract[wh,2],]~seq(2038,2068,by=1))
abline(a=lm(cout[extract[wh,1],extract[wh,2],]~seq(2038,2068,by=1))$coeff[1],b=lm(cout[extract[17,1],extract[17,2],]~seq(2038,2068,by=1))$coeff[2])

plot(hout[extract[wh,1],extract[wh,2],]~seq(2039,2068,by=1),ylab="HDD",xlab="Year",type="l",lwd=1.5,main="Trend in New York HDDs",
cex.lab=1.5,cex.main=1.5,cex.axis=1.5,cex=1.3)
points(hout[extract[wh,1],extract[wh,2],]~seq(2039,2068,by=1),pch=19)
abline(a=lm(hout[extract[wh,1],extract[wh,2],]~seq(2039,2068,by=1))$coeff[1],b=lm(hout[extract[17,1],extract[17,2],]~seq(2039,2068,by=1))$coeff[2])
points(hdd[extract[wh,1],extract[wh,2],]~seq(2039,2068,by=1))
dev.off()

## Here is some additional exploratory analysis stuff
## Looks like there is a coastal effect
seq1=seq(1,31,by=1)
i=1
plot(cdd[extract[i,1],extract[i,2],]~seq1, type="l",ylim=c(0,5000),col="white")
for (i in 1:30) {
  if(cities[i,6]=="M") {lines(cdd[extract[i,1],extract[i,2],]~seq1)} else {}
  if(cities[i,6]=="M") {text(seq1,cdd[extract[i,1],extract[i,2],],cities[i,6],)} else {}
}

seq2=seq(1,31,by=1)
i=1
plot(cdd[extract[i,1],extract[i,2],]~seq2, type="l",ylim=c(0,6000))
for (i in 2:30) {
  lines(cdd[extract[i,1],extract[i,2],]~seq2,col=i)
}
par(mfrow=c(2,2))
quilt.plot(loc,c(hdd[,1:85,1]),zlim=c(0,10000))
quilt.plot(loc,c(hdd[,1:85,11]),zlim=c(0,10000))
quilt.plot(loc,c(hdd[,1:85,21]),zlim=c(0,10000))
quilt.plot(loc,c(hdd[,1:85,30]),zlim=c(0,10000))

slope.comb=matrix(0,140,85)
se.slope.comb=matrix(0,140,85)
seq1=seq(1,30,by=1)
comb=hout+cout[,,1:30]

for (i in 1:140) {
  for (j in 1:85) {
    print(i+j/1000)
    slope.comb[i,j]=summary(lm(comb[i,j,]~seq1))$coeff[2,1]
    se.slope.comb[i,j]=summary(lm(comb[i,j,]~seq1))$coeff[2,2]
  }
}

cy2=matrix(0,30,4)
colnames(cy2)=c("Combined Trend","S.E.","Lower 95% Bound","Upper 95% Bound")
for (t in 1:30) {
  cy2[t,1]=round(slope.comb[extract[t,1],extract[t,2]],2)
  cy2[t,2]=round(se.slope.comb[extract[t,1],extract[t,2]],2)
  cy2[t,3]=round(slope.comb[extract[t,1],extract[t,2]]-conf*se.slope.cdd[extract[t,1],extract[t,2]],2)
  cy2[t,4]=round(slope.comb[extract[t,1],extract[t,2]]+conf*se.slope.cdd[extract[t,1],extract[t,2]],2)
}
data.frame(cities[,1],cy2)

par(mfrow=c(1,2))
quilt.plot(loc,c(slope.comb),nx=100,ny=85,main="Estimated Combined Trends")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
points(cities[,2],cities[,3],pch=19)

quilt.plot(loc,c(se.slope.comb),nx=100,ny=85,main="Estimated Combined Trends")
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
points(cities[,2],cities[,3],pch=19)

par(mfrow=c(1,2))
quilt.plot(loc,c(slope.comb)-2*c(se.slope.comb),nx=100,ny=85,main="Lower 95% Surface",zlim=c(-40,40))
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
points(cities[,2],cities[,3],pch=19)

quilt.plot(loc,c(slope.comb)+2*c(se.slope.comb),nx=100,ny=85,main="Upper 95% Surface",zlim=c(-40,40))
map(database="world", regions=c("usa","canada"), xlim=c(-150,-40), ylim=c(22, 60), lwd=2, add=T)
points(cities[,2],cities[,3],pch=19)




