##########################################
# FILE: EVModelingDJFSpatialVari0g.txt - WINTER
# AUTHOR: Elizabeth Shamseldin
# DATE: 4/13/06
# PURPOSE: Examine spatial trends in linear, quadratic, and cubic models (in lon, lat)
# Program CLEANS data for |xi|>1, logpsi, mu, grid not matching point, non-convergence
# 
#
# SPATIAL: Program to investigate spatial relationship.  
#			Looks at LogPt~Grid+Elev vs  LogPt~Grid+Elev+Lat+Lon vs  LogPt~Grid+Elev+Lat^2+Lon^2 vs LogPt~Grid+Elev+lat^3+lon^3
#			Investigates variograms for different models to determine if and when spatial trend is accounted for
#
# NOTE: On Elizabeth's Mac, X11 has to be RUNNING in order to produce png plots
#
# MODIFIED: 4/13/06
#
#		"NEW" param estimates downloaded from http://www.stat.unc.edu/postscript/rs/ext1/file3/
#		"NEW" estimates in EVdata
##########################################

StartTime <- Sys.time()

library(fields)
library(geoR)
library(evd)
library(ismev)

setwd("/Users/shamseld/Desktop/NCAR/")



#Set Confidence level z-value for plots
CI.z <- 1.96
#Set Header for Model Output
param.head <- cbind("1 intercept", "X1",  "aic",  "deviance", "numObs")
param.headE <- cbind("1 intercept", "X1",  "X2(Elev)", "aic",  "deviance", "numObs")
param.headELATLON <- cbind("1 intercept", "X1",  "X2(Elev)", "Lat", "Lon", "aic",  "deviance", "numObs")
param.headELATLONX <- cbind("1 intercept", "X1",  "X2(Elev)", "Lat", "Lon", "LatXLon", "aic",  "deviance", "numObs")


#Header
head <- cbind("num_station", "num_obs", "thresh", "lat", "lon", "elev", "conv", "mu_mle", "logpsi_mle", "xi_mle", "mu_se", "logpsi_se", "xi_se", "mulogpsi_corr", "muxi_corr", "logpsixi_corr")

#Read in grid data.  sep=""
g.DATA <- read.table(file="EVdata/GDJF19501999109511.txt", sep="", header=FALSE,  col.names=head) 
p.DATA <- read.table(file="EVdata/PDJF19501999109511.txt", sep="", header=FALSE,  col.names=head) 

StartTime <- Sys.time()
# Length: 4342
g.gridnum <- g.DATA[,1]
g.lat <- g.DATA[,4]
g.latmin <- g.lat-1.25
g.latmax <- g.lat+1.25
g.lon <- g.DATA[,5]
g.lonmin <- g.lon-1.25
g.lonmax <- g.lon+1.25
g.mu <- g.DATA[,8]
g.logpsi <- g.DATA[,9]
g.xi <- g.DATA[,10]
g.conv <- g.DATA[,7]

g.mu_se <- g.DATA[,11]
g.logpsi_se <- g.DATA[,12]
g.xi_se <- g.DATA[,13]
g.mulogpsi_corr <- g.DATA[,14]
g.muxi_corr <- g.DATA[,15]
g.logpsixi_corr <- g.DATA[,16]


p.lat <- p.DATA[,4]
p.lon <- p.DATA[,5]
p.elev <- p.DATA[,6]
p.mu <- p.DATA[,8]
p.logpsi <- p.DATA[,9]
p.xi <- p.DATA[,10]
p.gridnum <- matrix(0,length(p.DATA[,1]),1)
p.grid.lat <- matrix(0,length(p.DATA[,1]),1)
p.grid.lon <- matrix(0,length(p.DATA[,1]),1)

p.grid.mu <- matrix(0,length(p.DATA[,1]),1)
p.grid.logpsi <- matrix(0,length(p.DATA[,1]),1)
p.grid.xi <- matrix(0,length(p.DATA[,1]),1)
p.conv <- p.DATA[,7]
p.grid.conv <- matrix(0, length(p.DATA[,1]), 1)

p.grid.mu_se <- matrix(0,length(p.DATA[,1]),1)
p.grid.logpsi_se <- matrix(0,length(p.DATA[,1]),1)
p.grid.xi_se <- matrix(0,length(p.DATA[,1]),1)
p.grid.mulogpsi_corr <- matrix(0,length(p.DATA[,1]),1)
p.grid.muxi_corr <- matrix(0,length(p.DATA[,1]),1)
p.grid.logpsixi_corr <- matrix(0,length(p.DATA[,1]),1)

#Assign grids to point locations via lat and lon
for(i in 1:length(g.DATA[,1]))
	{
	for(j in 1:length(p.DATA[,1]))
		{
			if(p.lat[j] >= g.latmin[i] && p.lat[j] <= g.latmax[i] && p.lon[j] >= g.lonmin[i] && p.lon[j] <= g.lonmax[i])
			{p.gridnum[j] <- g.gridnum[i]}

		}#end j.p for
	}#end i.g for


#Assign mle values to point matrix
for(i in 1:length(p.DATA[,1]))
	{
		p.grid.lat[i] <- g.lat[(p.gridnum[i])]
		p.grid.lon[i] <- g.lon[(p.gridnum[i])]
		p.grid.mu[i] <- g.mu[(p.gridnum[i])]
		p.grid.logpsi[i] <- g.logpsi[(p.gridnum[i])]
		p.grid.xi[i] <- g.xi[(p.gridnum[i])]
		p.grid.conv[i] <- g.conv[(p.gridnum[i])]
		p.grid.mu_se[i] <- g.mu_se[(p.gridnum[i])]
		p.grid.logpsi_se[i] <- g.logpsi_se[(p.gridnum[i])]
		p.grid.xi_se[i] <- g.xi_se[(p.gridnum[i])]
		p.grid.mulogpsi_corr[i] <- g.mulogpsi_corr[(p.gridnum[i])]
		p.grid.muxi_corr[i] <- g.muxi_corr[(p.gridnum[i])]
		p.grid.logpsixi_corr[i] <- g.logpsixi_corr[(p.gridnum[i])]

	}


########################################################################################
########################################################################################
########################################################################################
#CLEANING DATA
# Program CLEANS data for |xi|>1, logpsi, mu, grid not matching point, non-convergence

Sys.time()  #Takes about 15 minutes

l <- length(p.DATA[,1])

p.TEMP <- cbind(p.DATA[1:l,1:16], p.gridnum[1:l], p.grid.conv[1:l], p.grid.lat[1:l], p.grid.lon[1:l], p.grid.mu[1:l], 
		p.grid.logpsi[1:l], p.grid.xi[1:l], p.grid.mu_se[1:l], p.grid.logpsi_se[1:l], p.grid.xi_se[1:l], 
		p.grid.mulogpsi_corr[1:l], p.grid.muxi_corr[1:l], p.grid.logpsixi_corr[1:l])


p.TEMP2 <- NULL

#for(i in 1:length(p.DATA[,1])){
for(i in 1:l){

	if( p.grid.xi[i] < 1 && p.grid.xi[i] > -1 && p.xi[i] < 1 && p.xi[i]  > -1 &&
		p.grid.logpsi[i] > 2 && p.grid.logpsi[i] < 6 && p.logpsi[i] > 2 && p.logpsi[i] < 6 &&
		p.grid.conv[i] == 0 && p.conv[i] == 0 ){

		p.TEMP2 <- rbind(p.TEMP2, p.TEMP[i,])
	}#end if

}#end for

Sys.time()

#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(p.TEMP2, file="PointVSGrid/WINTER/PDJFCl.csv")

########################################################################################
########################################################################################
########################################################################################
# Length: 4339
pc.lat <- p.TEMP2[,4]
pc.lon <- p.TEMP2[,5]
pc.elev <- p.TEMP2[,6]
pc.mu <- p.TEMP2[,8]
pc.logpsi <- p.TEMP2[,9]
pc.xi <- p.TEMP2[,10]
pc.mu_se <- p.TEMP2[,11]
pc.logpsi_se <- p.TEMP2[,12]
pc.xi_se <- p.TEMP2[,13]
pc.mulogpsi_corr <- p.TEMP2[,14]
pc.muxi_corr <- p.TEMP2[,15]
pc.logpsixi_corr <- p.TEMP2[,16]

pc.gridnum <- p.TEMP2[,17]
pc.grid.lat <- p.TEMP2[,19]
pc.grid.lon <- p.TEMP2[,20]

pc.grid.mu <- p.TEMP2[,21]
pc.grid.logpsi <- p.TEMP2[,22]
pc.grid.xi <- p.TEMP2[,23]
pc.conv <- p.TEMP2[,7]
pc.grid.conv <- p.TEMP2[,18]

pc.grid.mu_se <- p.TEMP2[,24]
pc.grid.logpsi_se <- p.TEMP2[,25]
pc.grid.xi_se <- p.TEMP2[,26]
pc.grid.mulogpsi_corr <- p.TEMP2[,27]
pc.grid.muxi_corr <- p.TEMP2[,28]
pc.grid.logpsixi_corr <- p.TEMP2[,29]


pc.return100 <- matrix(0, length(p.TEMP2[,1]), 1)
pc.grid.return100 <- matrix(0, length(p.TEMP2[,1]), 1)
pc.return100_se <- matrix(0, length(p.TEMP2[,1]), 1)
n.return <- 100

#Set return level based on xi=0, xi!=0
for(i in 1:length(p.TEMP2[,1]))
	{
		if(pc.xi[i] != 0) 
			{ pc.return100[i] <- (pc.mu[i] + exp(pc.logpsi[i])*((n.return^pc.xi[i] - 1)/pc.xi[i])) }
		if(pc.grid.xi[i] != 0)
 			{ pc.grid.return100[i] <- (pc.grid.mu[i] + exp(pc.grid.logpsi[i])*((n.return^pc.grid.xi[i] - 1)/pc.grid.xi[i])) }
		
		if(pc.xi[i] == 0)
			{ pc.return100[i] <- (pc.mu[i] + exp(pc.logpsi[i])*log(n.return)) }
 		if(pc.grid.xi[i] == 0)
 			{ pc.grid.return100[i] <- (pc.grid.mu[i] + exp(pc.grid.logpsi[i])*log(n.return)) }
 		
 		pc.return100_se[i] <- delta (pc.mu[i], pc.logpsi[i], pc.xi[i], pc.mu_se[i], pc.logpsi_se[i], pc.xi_se[i],
					 pc.mulogpsi_corr[i], pc.muxi_corr[i], pc.logpsixi_corr[i], n.return) 
	}


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


return100LOGelevC.glm <- glm(log(pc.return100) ~ pc.grid.return100 + pc.elev) 

########################################################################
#Clean studentized residuals magnitude greater than 2
pc.return.weightE <- matrix(1, length(p.TEMP2[,1]), 1)

for(i in 1:length(p.TEMP2[,1]))  #This takes ~10 minutes to run
	{
		if(rstudent(return100LOGelevC.glm)[i] > 2 || rstudent(return100LOGelevC.glm)[i] < -2)
		{pc.return.weightE[i] <- 0}

	}

pc.grid.return100WE <- NULL
pc.return100WE <- NULL
pc.elev100WE <- NULL
pc.lat100WE <- NULL
pc.lon100WE <- NULL
pc.return100_seW <- NULL


for(i in 1:length(p.TEMP2[,1]))  
	{
		if(pc.return.weightE[i] != 0)
		{pc.grid.return100WE <- rbind(pc.grid.return100WE, pc.grid.return100[i]);
		pc.return100WE <- rbind(pc.return100WE, pc.return100[i]);
		pc.elev100WE <- rbind(pc.elev100WE, pc.elev[i]);
		pc.lat100WE <- rbind(pc.lat100WE, pc.lat[i]);
		pc.lon100WE <- rbind(pc.lon100WE, pc.lon[i]);
		pc.return100_seW <- rbind(pc.return100_seW, pc.return100_se[i])	}

	}
#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
#write.table(cbind(pc.grid.return100WE, pc.return100WE, pc.elevWE), file="PointVSGrid/WINTER/PDJF100WE.csv")

Sys.time()
###################################################
###################################################
return100LOGelevCW.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE)  
pred100LOGElevW <- predict(return100LOGelevCW.glm, se.fit=TRUE)


########################################################################################
########################################################################################
########################################################################################
# VARIOGRAMS

#Set max.dist to 1000?
vgram.v <- vgram(cbind(pc.lon, pc.lat), return100LOGelevC.glm$residuals, lon.lat=TRUE, breaks=seq(0,3000,100))


png(filename="PointVSGrid/WINTER/Variogram/LogPt+ElevVgram-DJF95.png")
plot(vgram.v$center, vgram.v$stats[2,], main ="Vgram of LOGelev - DJF 95")
dev.off()

vgram.vMD1000 <- vgram(cbind(pc.lon, pc.lat), return100LOGelevC.glm$residuals, lon.lat=TRUE, breaks=seq(0,1000,50), dmax=1000)

png(filename="PointVSGrid/WINTER/Variogram/LogPt+ElevVgramMD1000-DJF95.png")
plot(c(0,1000), c(0,0.25), type="n", main ="Vgram of LOGelev MD1000 - DJF 95", xlab="Distance", ylab="vgm")
points(vgram.vMD1000$center, vgram.vMD1000$stats[2,])
dev.off()

###################################################
pc.lat100WE2 <- pc.lat100WE^2
pc.lon100WE2 <- pc.lon100WE^2
#AIC 1496 DJF
return100LOGelevLATLON2.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE
									+ pc.lat100WE*pc.lon100WE + pc.lat100WE2 + pc.lon100WE2)  
pred100LOGelevLATLON2 <- predict(return100LOGelevLATLON2.glm, se.fit=TRUE)

									
vgram.vLOGelevLATLON2 <- vgram(cbind(pc.lon, pc.lat), return100LOGelevLATLON2.glm$residuals, lon.lat=TRUE, breaks=seq(0,3000,100))


png(filename="PointVSGrid/WINTER/Variogram/LogPt+ElevLATLON2Vgram-DJF95.png")
plot(vgram.vLOGelevLATLON2$center, vgram.vLOGelevLATLON2$stats[2,], main ="Vgram of LOGelevLATLON2 - DJF 95")
dev.off()

vgram.vLOGelevLATLON2md1000 <- vgram(cbind(pc.lon, pc.lat), return100LOGelevLATLON2.glm$residuals, lon.lat=TRUE, 
	breaks=seq(0,1000,50), dmax=1000)

png(filename="PointVSGrid/WINTER/Variogram/LogPt+ElevLATLON2VgramMD1000-DJF95.png")
plot(c(0,1000), c(0,0.25), type="n", xlab="Distance", ylab="vgm", main ="Vgram of LOGelevLATLON2 MD1000 - DJF 95")
points(vgram.vLOGelevLATLON2md1000$center, vgram.vLOGelevLATLON2md1000$stats[2,])
dev.off()



png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevLATLON2StuResidDJF95.png")
plot(pc.grid.return100WE, rstudent(return100LOGelevLATLON2.glm), main="100-Yr Ret LogPt+Elev+LAT2+LON2 StuResid DJFW", cex=0.4)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElev+LAT2LON2FitDJFW.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Ret LogPt+Elev+LAT2+LON2 Fitted DJF95", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGelevLATLON2.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGelevLATLON2.glm$fitted.values[i]+CI.z*pred100LOGelevLATLON2$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGelevLATLON2.glm$fitted.values[i]-CI.z*pred100LOGelevLATLON2$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

###
#R^2
SSEquad <-  sum((return100LOGelevLATLON2.glm$fitted.values-log(pc.return100WE))*(return100LOGelevLATLON2.glm$fitted.values-log(pc.return100WE)))
SSY <-  sum((mean(log(pc.return100WE))-log(pc.return100WE))*(mean(log(pc.return100WE))-log(pc.return100WE)))
R2quad <- 1-SSEquad/SSY  #0.7440804


###################################################
pc.lat100WE3 <- pc.lat100WE^3
pc.lon100WE3 <- pc.lon100WE^3

#AIC  1062 DJF
return100LOGelevLATLON3.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE
									+ pc.lat100WE*pc.lon100WE
									+ pc.lat100WE2 + pc.lon100WE2 + pc.lon100WE*pc.lat100WE2 + pc.lat100WE*pc.lon100WE2 
									+ pc.lat100WE3 + pc.lon100WE3)  
pred100LOGelevLATLON3 <- predict(return100LOGelevLATLON3.glm, se.fit=TRUE)
# min(pred100LOGelevLATLON3$se.fit)
# 0.0076593
# max(pred100LOGelevLATLON3$se.fit)
# 0.054475


vgram.vLOGelevLATLON3 <- vgram(cbind(pc.lon, pc.lat), return100LOGelevLATLON3.glm$residuals, lon.lat=TRUE, breaks=seq(0,3000,100))

png(filename="PointVSGrid/WINTER/Variogram/LogPt+ElevLATLON3Vgram-DJF95.png")
plot(vgram.vLOGelevLATLON3$center, vgram.vLOGelevLATLON3$stats[2,], main ="Vgram of LOGelevLATLON3 - DJF 95")
dev.off()

vgram.vLOGelevLATLON3md1000 <- vgram(cbind(pc.lon, pc.lat), return100LOGelevLATLON3.glm$residuals, lon.lat=TRUE, 
	breaks=seq(0,1000,50), dmax=1000)
	
png(filename="PointVSGrid/WINTER/Variogram/LogPt+ElevLATLON3VgramMD1000-DJF95.png")
plot(c(0,1000), c(0,0.25), type="n", xlab="Distance", ylab="vgm", main ="Vgram of LOGelevLATLON3 MD1000 - DJF 95")
points(vgram.vLOGelevLATLON3md1000$center, vgram.vLOGelevLATLON3md1000$stats[2,])
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevLATLON3StuResidDJF95.png")
plot(pc.grid.return100WE, rstudent(return100LOGelevLATLON3.glm), main="100-Yr Ret LogPt+Elev+LAT3+LON3 StuResid DJFW", cex=0.4)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElev+LAT3LON3FitDJFW.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Ret LogPt+Elev+LAT3+LON3 Fitted DJF95", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGelevLATLON3.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGelevLATLON3.glm$fitted.values[i]+CI.z*pred100LOGelevLATLON3$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGelevLATLON3.glm$fitted.values[i]-CI.z*pred100LOGelevLATLON3$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

###
#R^2
SSEcub <-  sum((return100LOGelevLATLON3.glm$fitted.values-log(pc.return100WE))*(return100LOGelevLATLON3.glm$fitted.values-log(pc.return100WE)))
R2cub <- 1-SSEcub/SSY  #0.770668

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

png(filename="PointVSGrid/WINTER/Variogram/VariogCompare-DJF95.png")
plot(c(0,1000), c(0,0.25), type="n", xlab="Distance", ylab="vgm", main ="Vgram Comparison MD1000 - DJF 95")
points(vgram.vMD1000$center, vgram.vMD1000$stats[2,])
points(vgram.vLOGelevLATLON2md1000$center, vgram.vLOGelevLATLON2md1000$stats[2,], cex=0.7)
points(vgram.vLOGelevLATLON3md1000$center, vgram.vLOGelevLATLON3md1000$stats[2,], cex=0.4)
dev.off()


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

FITDIFFquadcub <- exp(return100LOGelevLATLON2.glm$fitted.values)-exp(return100LOGelevLATLON3.glm$fitted.values)
maxFITDIFF <- max(FITDIFFquadcub) #534
aveFITDIFF <- mean(FITDIFFquadcub) # -1.421194
minFITDIFF <- min(FITDIFFquadcub) # -447

loc <- cbind(pc.lon100WE, pc.lat100WE)

 aiFITDIFF <- as.image(FITDIFFquadcub, x = loc)
 
png(filename="PointVSGrid/WINTER/ImagePlots/FITDIFFquadcub-DJF95.png")
  image.plot(aiFITDIFF, xlab="longitude", ylab="latitude", main="Diff in Quad-Cubic Fitted (Ave=-1.421194) DJF-95")
  US(add=T)
  dev.off()
####
FITRATIOquadcub <- exp(return100LOGelevLATLON2.glm$fitted.values)/exp(return100LOGelevLATLON3.glm$fitted.values)

maxFITRATIO <- max(FITRATIOquadcub) #1.442546
aveFITRATIO <- mean(FITRATIOquadcub) # 1.004392
minFITRATIO <- min(FITRATIOquadcub) # 0.6769054

loc <- cbind(pc.lon100WE, pc.lat100WE)

 aiFITRATIO <- as.image(FITRATIOquadcub, x = loc)
 
png(filename="PointVSGrid/WINTER/ImagePlots/FITRATIOquadcub-DJF95.png")
  image.plot(aiFITRATIO, xlab="longitude", ylab="latitude", main="Ratio of Quad/Cubic Fitted (Ave=1.004392) DJF-95")
  US(add=T)
  dev.off()
  
########
loc <- cbind(pc.lon100WE, pc.lat100WE)


SElogElev <- exp(pred100LOGElevW$se.fit)

 aiSElogElev <- as.image(SElogElev, x = loc)
 
png(filename="PointVSGrid/WINTER/ImagePlots/SElogElev-DJF95.png")
  image.plot(aiSElogElev, xlab="longitude", ylab="latitude", main="SE of LogPt+Elev Fitted DJF-95")
  US(add=T)
  dev.off()
###
SEquad <- exp(pred100LOGelevLATLON2$se.fit)
 aiSEquad <- as.image(SEquad, x = loc)
 maxSEquad <- max(SEquad) #1.037356

png(filename="PointVSGrid/WINTER/ImagePlots/SEquad-DJF95.png")
  image.plot(aiSEquad, xlab="longitude", ylab="latitude", main="SE of Quad Fitted DJF-95")
  US(add=T)
  dev.off()
###  
SEcub <- exp(pred100LOGelevLATLON3$se.fit)  
aiSEcub <- as.image(SEcub, x = loc)
 
png(filename="PointVSGrid/WINTER/ImagePlots/SEcub-DJF95.png")
  image.plot(aiSEcub, xlab="longitude", ylab="latitude", main="SE of Cubic Fitted DJF-95")
  US(add=T)
  dev.off()
  
  
####################################################
####################################################
####################################################
#Investigating Standard Errors:

plot(pc.return100WE, pc.return100_seW, xlab="Point Return", ylab="Std Err", main="Std Error vs Point Return")

png(filename="PointVSGrid/WINTER/Variogram/StdErrVsFITDIFFquadcub-DJF95.png")
plot(pc.return100_seW, abs(FITDIFFquadcub), xlab="Pt StdErr", ylab="Abs(Fit Diff): Quad-Cubic", 
		main="Abs(Fit Diff) VS Pt-100Ret Std Err: Quad-Cubic - DJF95", cex=0.4)
		lines(pc.return100_seW,pc.return100_seW)
dev.off()

aiSEreturnLOG <- as.image(log(pc.return100_seW), x = loc)
png(filename="PointVSGrid/WINTER/ImagePlots/LOGSEreturn-DJF95.png")
  image.plot(aiSEreturnLOG, xlab="longitude", ylab="latitude", main="LOG of SE of Point 100-yr Return DJF-95")
  US(add=T)
  dev.off()
  
  aiSEreturn <- as.image(pc.return100_seW, x = loc)
png(filename="PointVSGrid/WINTER/ImagePlots/SEreturn-DJF95.png")
  image.plot(aiSEreturn, zlim=c(0,600), xlab="longitude", ylab="latitude", main="SE of Point 100-yr Return DJF-95")
  US(add=T)
  dev.off()
  

  
####################################################
####################################################
####################################################
#Quartic Model:

pc.lat100WE4 <- pc.lat100WE^4
pc.lon100WE4 <- pc.lon100WE^4

#AIC  614.4 DJF
return100LOGelevLATLON4.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE
									+ pc.lat100WE*pc.lon100WE
									+ pc.lat100WE2 + pc.lon100WE2 + pc.lon100WE*pc.lat100WE2 + pc.lat100WE*pc.lon100WE2 
									+ pc.lat100WE3 + pc.lon100WE3
									+ pc.lat100WE4 + pc.lat100WE3*pc.lon100WE + pc.lat100WE2*pc.lon100WE2
									+ pc.lat100WE*pc.lon100WE3 + pc.lon100WE4)  

pred100LOGelevLATLON4 <- predict(return100LOGelevLATLON4.glm, se.fit=TRUE)


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevLATLON4StuResidDJFW.png")
plot(pc.grid.return100WE, rstudent(return100LOGelevLATLON4.glm), main="100-Yr Ret LogPt+Elev+LAT4+LON4 StuResid DJF95", cex=0.4)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElev+LAT4LON4FitDJFW.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Ret LogPt+Elev+LAT4+LON4 Fitted DJF95", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGelevLATLON4.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGelevLATLON4.glm$fitted.values[i]+CI.z*pred100LOGelevLATLON4$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGelevLATLON4.glm$fitted.values[i]-CI.z*pred100LOGelevLATLON4$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

###
#R^2
SSEquart <-  sum((return100LOGelevLATLON4.glm$fitted.values-log(pc.return100WE))*(return100LOGelevLATLON4.glm$fitted.values-log(pc.return100WE)))
R2quart <- 1-SSEquart/SSY  #0.795327

####################################################
####################################################
####################################################
#NOTES ON AIC CALC IN R:
#This is a generic function, with methods in base R for "aov", "coxph", "glm", "lm", "negbin" and "survreg" classes.
#The criterion used is
#AIC = - 2*log L + k * edf,
# where L is the likelihood and edf the equivalent degrees of freedom 
#(i.e., the number of parameters for usual parametric models) of fit.

#For linear models with unknown scale (i.e., for lm and aov), -2log L is computed from the deviance 
#and uses a different additive constant to AIC.
#k = 2 corresponds to the traditional AIC, using k =  log(n) provides the BIC (Bayes IC) instead.

#http://www.r-project.org/nocvs/mail/r-help/2001/3337.html
#AIC = n ln(sigma^2) + 2*parms, 
# where sigma^2 is the estimate for MSE? 

#http://economics.about.com/cs/economicsglossary/g/akaikes.htm
#AIC=T ln(RSS) + 2K 
#www.stat.wisc.edu/course/ st333-larget/public/html/aic.pdf
#AIC = n + n*log(2pi) + n*log(MSE) +2*(p+1)
#BIC = n + n*log(2pi) + n*log(MSE) +log(n)*(p+1)

#########
#CALCULATION OF AIC USING n*ln(MSE)+2*k

#log(Pt) + Elev #AIC R: 3424.572
diff <- return100LOGelevCW.glm$fitted.values-log(pc.return100WE)
diffSQ <- diff*diff
AICelev <- length(return100LOGelevCW.glm$fitted.values)*log(mean(diffSQ)) + 2*length(return100LOGelevCW.glm$coefficients)

#AICelev = -7999.883



#log(Pt) + Elev + LAT2 + LON2 + X #AIC R: 1495.772
diffLATLON2 <- return100LOGelevLATLON2.glm$fitted.values-log(pc.return100WE)
diffLATLON2SQ <- diffLATLON2*diffLATLON2
AICLATLON2 <- length(return100LOGelevLATLON2.glm$fitted.values)*log(mean(diffLATLON2SQ))+2*length(return100LOGelevLATLON2.glm$coefficients)
#AICLATLON2 = -9928.683
BICLATLON2 <- length(return100LOGelevLATLON2.glm$fitted.values)*log(mean(diffLATLON2SQ))+
	length(return100LOGelevLATLON2.glm$fitted.values)*(log(length(return100LOGelevLATLON2.glm$coefficients))+1)
#BICLATLON2 = 2450.069
#2450.069 + 4025 + 4025*log(2*pi) = 13872.52


#log(Pt) + Elev + LAT3 + LON3 + X #AIC R: 1062.260
diffLATLON3 <- return100LOGelevLATLON3.glm$fitted.values-log(pc.return100WE)
diffLATLON3SQ <- diffLATLON3*diffLATLON3
AICLATLON3 <- length(return100LOGelevLATLON3.glm$fitted.values)*log(mean(diffLATLON3SQ))+2*length(return100LOGelevLATLON3.glm$coefficients)
#AICLATLON3 = -10362.19
# -10362.19 + 2 + 4025 + 4025*log(2*pi) = 1062.265
BICLATLON3 <- length(return100LOGelevLATLON3.glm$fitted.values)*log(mean(diffLATLON3SQ))+
		length(return100LOGelevLATLON3.glm$fitted.values)*(log(length(return100LOGelevLATLON3.glm$coefficients))+1)
#BICLATLON3 = 3640.555
#3640.555 + 4025 + 4025*log(2*pi) = 15063.01

#log(Pt) + Elev + LAT4 + LON4 + X #AIC R: 614.3872
diffLATLON4 <- return100LOGelevLATLON4.glm$fitted.values-log(pc.return100WE)
diffLATLON4SQ <- diffLATLON4*diffLATLON4
AICLATLON4 <- length(return100LOGelevLATLON4.glm$fitted.values)*log(mean(diffLATLON4SQ))+2*length(return100LOGelevLATLON4.glm$coefficients)
#AICLATLON4 = -10810.07
# -10810.07 + 2 + 4025 + 4025*log(2*pi) = 614.3852
BICLATLON4 <- length(return100LOGelevLATLON4.glm$fitted.values)*log(mean(diffLATLON4SQ))+
			length(return100LOGelevLATLON4.glm$fitted.values)*(log(length(return100LOGelevLATLON4.glm$coefficients))+1)
#BICLATLON4 = 4584.616
#4584.616 + 4025 + 4025*log(2*pi) = 16007.07

FcubVSquartic <- ((sum(diffLATLON3SQ) - sum(diffLATLON4SQ))*return100LOGelevLATLON4.glm$df.residual)/
					((return100LOGelevLATLON3.glm$df.residual - return100LOGelevLATLON4.glm$df.residual)*sum(diffLATLON4SQ))
#FcubVSquartic = 96.57693