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


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


return100LOG.nls <- nls(log(pc.return100) ~ b0+b1*pc.grid.return100^b2, start=list(b0=0,b1=1,b2=1)) 

pred100LOGnls <- predict(return100LOG.nls, se.fit=TRUE)
#out100LOGnls <- c()
#write.table(t(out100LOGnls), file="PointVSGrid/WINTER/Models/NLS/100Ret/ModelParams/100RetDJFLOGnls.csv", 
#col.names=param.headNLS)



png(filename="PointVSGrid/WINTER/Models/NLS/100Ret/PNG/100RetLogNLSResidDJFW.png")
plot(pc.grid.return100, (pc.return100-pred100LOGnls), main="100-Yr Ret LogPt NLS Resid DJF", cex=0.4)
dev.off()

#sdPred100LOGnls <- sd(pred100LOGnls)
#png(filename="PointVSGrid/WINTER/Models/NLS/100Ret/PNG/100RetLogNLSStandardizedResidDJFW.png")
#plot(pc.grid.return100, (pc.return100-pred100LOGnls)/sdPred100LOGnls, main="100-Yr Ret LogPt NLS StandardizedResid DJF", cex=0.4)
#dev.off()
 
 
 loc <- cbind(pc.lon, pc.lat)
 aiFitNLS <- as.image(pred100LOGnls, x = loc)

png(filename="PointVSGrid/WINTER/ImagePlots/LogPtNLSFitted-DJF95.png")
  image.plot(aiFitNLS, xlab="longitude", ylab="latitude", main="LogPt~Grid NLS FItted DJF-95")
  US(add=T)
  dev.off()

 aiResid <- as.image((pc.return100-pred100LOGnls), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPtNLSResid-DJF95.png")
  image.plot(aiResid, zlim=c(-200,2000), xlab="longitude", ylab="latitude", main="LogPt~Grid NLS Resid DJF-95")
  US(add=T)
  dev.off()

 


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

return100LOGElev.nls <- nls(log(pc.return100) ~ b0+b1*pc.grid.return100^b2+b3*pc.elev, start=list(b0=0,b1=1,b2=1,b3=-0.0001)) 

pred100LOGElevnls <- predict(return100LOGElev.nls, se.fit=TRUE)
#out100LOGnls <- c()
#write.table(t(out100LOGnls), file="PointVSGrid/WINTER/Models/NLS/100Ret/ModelParams/100RetDJFLOGnls.csv", 
#col.names=param.headNLS)


png(filename="PointVSGrid/WINTER/Models/NLS/100Ret/PNG/100RetLogElevNLSResidDJFW.png")
plot(pc.grid.return100, (pc.return100-pred100LOGElevnls), main="100-Yr Ret LogPt+Elev NLS Resid DJF", cex=0.4)
dev.off()

 
 aiFitENLS <- as.image(pred100LOGElevnls, x = loc)

png(filename="PointVSGrid/WINTER/ImagePlots/LogPtElevNLSFitted-DJF95.png")
  image.plot(aiFitENLS, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev NLS FItted DJF-95")
  US(add=T)
  dev.off()



 aiResidE <- as.image((pc.return100-pred100LOGElevnls), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPtElevNLSResid-DJF95.png")
  image.plot(aiResidE,  zlim=c(-200,2000), xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev NLS Resid DJF-95")
  US(add=T)
  dev.off()

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

return100LOGElevLATLON.nls <- nls(log(pc.return100) ~ b0+b1*pc.grid.return100^b2+b3*pc.elev+b4*pc.lat+b5*pc.lon, 
								start=list(b0=0,b1=1,b2=1,b3=-0.0001, b4=1,b5=1)) 

pred100LOGElevLATLONnls <- predict(return100LOGElevLATLON.nls, se.fit=TRUE)


png(filename="PointVSGrid/WINTER/Models/NLS/100Ret/PNG/100RetLogElevLATLONNLSResidDJFW.png")
plot(pc.grid.return100, (pc.return100-pred100LOGElevLATLONnls), main="100-Yr Ret LogPt+Elev+Lat+Lon NLS Resid DJF", cex=0.4)
dev.off()

 
 aiFitELATLONnls <- as.image(pred100LOGElevLATLONnls, x = loc)

png(filename="PointVSGrid/WINTER/ImagePlots/LogPtElevLATLONNLSFitted-DJF95.png")
  image.plot(aiFitELATLONnls, xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon NLS FItted DJF-95")
  US(add=T)
  dev.off()



 aiResidELATLON <- as.image((pc.return100-pred100LOGElevLATLONnls), x = loc)
 
  png(filename="PointVSGrid/WINTER/ImagePlots/LogPtElevLATLONNLSResid-DJF95.png")
  image.plot(aiResidELATLON, zlim=c(-200,2000), xlab="longitude", ylab="latitude", main="LogPt~Grid+Elev+Lat+Lon NLS Resid DJF-95")
  US(add=T)
  dev.off()