##########################################
# FILE: EVModelingDJFSpatial.txt - WINTER
# AUTHOR: Elizabeth Shamseldin
# DATE: 3/27/06
# PURPOSE: Create point-level yr-return estimates modeled from grid-level yr-returns
# using parameter estimates from
# http://www.stat.unc.edu/postscript/rs/ext1/ for NCAR point vs grid-cell models
# 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 resid pattern.  Expect to see spatial pattern across resids.
#			Looks at n.return spatial structure across grid.  
#			Investigates effect of including LAT, LON in model.
#			Looks at that residual pattern - hopefully no spatial trend is seen, lat-lon accounted for it.
#
# NOTE: On Elizabeth's Mac, X11 has to be RUNNING in order to produce png plots
#
# MODIFIED: 3/27/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(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.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)
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)) }
 		
	}


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


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


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])	}

	}
#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)  


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

aiReturnPT <- as.image(pc.return100WE, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/Point100Return-DJF95.png")
  image.plot(aiReturnPT, zlim=c(100,1500), xlab="longitude", ylab="latitude", main="Point 100-Yr Return DJF-95")
  US(add=T)
  dev.off()
  
  aiReturnGRID <- as.image(pc.grid.return100WE, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/Grid100Return-DJF95.png")
  image.plot(aiReturnGRID, xlab="longitude", ylab="latitude", main="Grid 100-Yr Return DJF-95")
  US(add=T)
  dev.off()

 aiStuResid <- as.image(rstudent(return100LOGelevCW.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevStuResid-DJF95.png")
  image.plot(aiStuResid, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev StuResid DJF-95")
  US(add=T)
  dev.off()

 aiFitted <- as.image(return100LOGelevCW.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevFitted-DJF95.png")
  image.plot(aiFitted, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev Fitted DJF-95")
  US(add=T)
  dev.off()

#AIC 1922 DJF
return100LOGelevLATLON.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE)  

pred100LOGElevLATLON <- predict(return100LOGelevLATLON.glm, se.fit=TRUE)
out100LOGElevLATLON <- c(return100LOGelevLATLON.glm$coefficients, return100LOGelevLATLON.glm$aic, return100LOGelevLATLON.glm$deviance,
length(return100LOGelevLATLON.glm$residuals))
write.table(t(out100LOGElevLATLON), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtElevModelLATLON.csv", 
col.names=param.headELATLON)

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


 aiFittedLATLON <- as.image(return100LOGelevLATLON.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+Elev+Lat+LonFittedDJF95.png")
  image.plot(aiFittedLATLON, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon Fitted DJF-95")
  US(add=T)
  dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevLatLonFitDJFW.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Ret LogPt+Elev+Lat+Lon Fitted Values DJF-W", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGelevLATLON.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGelevLATLON.glm$fitted.values[i]+CI.z*pred100LOGElevLATLON$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGelevLATLON.glm$fitted.values[i]-CI.z*pred100LOGElevLATLON$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()


 aiStuResidLATLON <- as.image(rstudent(return100LOGelevLATLON.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLonStuResid-DJF95.png")
  image.plot(aiStuResidLATLON, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon StuResid DJF-95")
  US(add=T)
  dev.off()




#AIC 1703 DJF
return100LOGelevLATLONX.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE
									+ pc.lat100WE*pc.lon100WE)  
pred100LOGElevLATLONX <- predict(return100LOGelevLATLONX.glm, se.fit=TRUE)
out100LOGElevLATLONX <- c(return100LOGelevLATLONX.glm$coefficients, return100LOGelevLATLONX.glm$aic, return100LOGelevLATLONX.glm$deviance,
length(return100LOGelevLATLONX.glm$residuals))
write.table(t(out100LOGElevLATLONX), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtElevModelLATLONX.csv", 
col.names=param.headELATLONX)


aiStuResidLATLONX <- as.image(rstudent(return100LOGelevLATLONX.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLonXStuResid-DJF95.png")
  image.plot(aiStuResidLATLONX, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon+X StuResid DJF-95")
  US(add=T)
  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)
out100LOGElevLATLON2 <- c(return100LOGelevLATLON2.glm$coefficients, return100LOGelevLATLON2.glm$aic, return100LOGelevLATLON2.glm$deviance,
length(return100LOGelevLATLON2.glm$residuals))
#write.table(t(out100LOGElevLATLON2), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtElevModelLATLON2.csv", 
col.names=param.headELATLON)



 aiFittedLATLON2 <- as.image(return100LOGelevLATLON2.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+Elev+Lat2+Lon2FittedDJF95.png")
  image.plot(aiFittedLATLON2, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat2+Lon2 Fitted DJF-95")
  US(add=T)
  dev.off()
  
aiStuResidLATLON2 <- as.image(rstudent(return100LOGelevLATLON2.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLon2StuResid-DJF95.png")
  image.plot(aiStuResidLATLON2, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat2+Lon2 StuResid DJF-95")
  US(add=T)
  dev.off()
  
  

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)
out100LOGElevLATLON3 <- c(return100LOGelevLATLON3.glm$coefficients, return100LOGelevLATLON3.glm$aic, return100LOGelevLATLON3.glm$deviance,
length(return100LOGelevLATLON3.glm$residuals))
#write.table(t(out100LOGElevLATLON3), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtElevModelLATLON3.csv", 
col.names=param.headELATLON)



 aiFittedLATLON3 <- as.image(return100LOGelevLATLON3.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+Elev+Lat3+Lon3FittedDJF95.png")
  image.plot(aiFittedLATLON3, zlim=c(5.3,7.7), xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat3+Lon3 Fitted DJF-95")
  US(add=T)
  dev.off()

aiStuResidLATLON3 <- as.image(rstudent(return100LOGelevLATLON3.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLon3StuResid-DJF95.png")
  image.plot(aiStuResidLATLON3, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat3+Lon3 StuResid 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)


 aiFittedLATLON4 <- as.image(return100LOGelevLATLON4.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+Elev+Lat4+Lon4FittedDJF95.png")
  image.plot(aiFittedLATLON4, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat4+Lon4 Fitted DJF-95")
  US(add=T)
  dev.off()

aiStuResidLATLON4 <- as.image(rstudent(return100LOGelevLATLON4.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLon4StuResid-DJF95.png")
  image.plot(aiStuResidLATLON4, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat4+Lon4 StuResid DJF-95")
  US(add=T)
  dev.off()







#############################
#Quadratic in Elevation

pc.elev100WE2 <- pc.elev100WE^2

#AIC 1679 DJF
return100LOGelev2LATLON.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE
									+ pc.elev100WE2)  

pred100LOGElev2LATLON <- predict(return100LOGelev2LATLON.glm, se.fit=TRUE)
out100LOGElev2LATLON <- c(return100LOGelev2LATLON.glm$coefficients, return100LOGelev2LATLON.glm$aic, 
return100LOGelev2LATLON.glm$deviance, length(return100LOGelev2LATLON.glm$residuals))


 aiFittedLATLONElev2 <- as.image(return100LOGelev2LATLON.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+Elev+Lat+Lon+Elev2FittedDJF95.png")
  image.plot(aiFittedLATLONElev2, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon+Elev2 Fitted DJF-95")
  US(add=T)
  dev.off()
  
aiStuResidLATLONElev2 <- as.image(rstudent(return100LOGelev2LATLON.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLonElev2StuResid-DJF95.png")
  image.plot(aiStuResidLATLONElev2, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon+Elev2 StuResid DJF-95")
  US(add=T)
  dev.off()
  
#############
#AIC  DJF
return100LOGelev2LATLONX.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elev100WE + pc.lat100WE + pc.lon100WE
									+ pc.elev100WE2 + pc.elev100WE*pc.lat100WE + pc.elev100WE*pc.lon100WE)  

pred100LOGElev2LATLONX <- predict(return100LOGelev2LATLONX.glm, se.fit=TRUE)
out100LOGElev2LATLONX <- c(return100LOGelev2LATLONX.glm$coefficients, return100LOGelev2LATLONX.glm$aic, 
return100LOGelev2LATLONX.glm$deviance, length(return100LOGelev2LATLONX.glm$residuals))


 aiFittedLATLONElev2X <- as.image(return100LOGelev2LATLONX.glm$fitted.values, x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+Elev+Lat+Lon+Elev2XFittedDJF95.png")
  image.plot(aiFittedLATLONElev2X, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon+Elev2+X Fitted DJF-95")
  US(add=T)
  dev.off()
  
aiStuResidLATLONElev2X <- as.image(rstudent(return100LOGelev2LATLONX.glm), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPt+ElevLatLonElev2XStuResid-DJF95.png")
  image.plot(aiStuResidLATLONElev2X, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon+Elev2+X StuResid DJF-95")
  US(add=T)
  dev.off()
  



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

pc.return50 <- matrix(0, length(p.TEMP2[,1]), 1)
pc.grid.return50 <- matrix(0, length(p.TEMP2[,1]), 1)
n.return50 <- 50

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




return50LOGelevC.glm <- glm(log(pc.return50) ~ pc.grid.return50 + pc.elev)