##########################################
# FILE: ImagePlotDFJF97.txt
# AUTHOR: Elizabeth Shamseldin
# DATE: 3/21/06
# PURPOSE: Create image plots of parameter estimates for comparison
#	across grid vs point and 95 vs 97 thresholds, seasonally
#
# MODIFIED: 
##########################################

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

setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER/ImagePlots/97")


#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="/Users/shamseld/Desktop/NCAR/EVdata/GDJF19701999109711.txt", sep="", header=FALSE,  col.names=head) #WINTERgrid
p.DATA <- read.table(file="/Users/shamseld/Desktop/NCAR/EVdata/PDJF19701999109711.txt", sep="", header=FALSE,  col.names=head) #WINTERgrid



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.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.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)

### Grid IMAGE.PLOTS
 
 mtxMU <- matrix(g.mu, ncol = 24)
  mtxMU <- t(mtxMU)
 png(filename="GridMuDJF97.png")
  image.plot(mtxMU, x = unique(g.lon), y = unique(g.lat), xlab="longitude", ylab="latitude", main="Grid Mu DJF-97")
  US(add = T)
 dev.off()
 
 mtxLOGPSI <- matrix(g.logpsi, ncol = 24)
 mtxLOGPSI <- t(mtxLOGPSI)
 png(filename="GridLogPsiDJF97.png")
  image.plot(mtxLOGPSI, x = unique(g.lon), y = unique(g.lat), xlab="longitude", ylab="latitude", main="Grid LogPsi DJF-97")
  US(add = T)
 dev.off()

 
  mtxXI <- matrix(g.xi, ncol = 24)
  mtxXI <- t(mtxXI)
 png(filename="GridXiDJF97.png")
  image.plot(mtxXI, x = unique(g.lon), y = unique(g.lat), xlab="longitude", ylab="latitude", main="Grid Xi DJF-97")
  US(add = T)
 dev.off()
  

 

### Point IMAGE.PLOTS

 loc <- cbind(p.lon, p.lat)

 aiMU <- as.image(p.mu, x = loc)
 
  png(filename="PointMuDJF97.png")
  image.plot(aiMU, zlim=c(0,800), xlab="longitude", ylab="latitude", main="Point Mu DJF-97")
  US(add=T)
  dev.off()

 aiLOGPSI <- as.image(p.logpsi, x = loc)

 png(filename="PointLogPsiDJF97.png")
   image.plot(aiLOGPSI, zlim=c(3,5), xlab="longitude", ylab="latitude", main="Point LogPsi DJF-97" )
  US(add=T)
  dev.off()
  
  aiXI <- as.image(p.xi, x = loc)

 png(filename="PointXiDJF97.png")
   image.plot(aiXI, zlim=c(-0.2,0.8), xlab="longitude", ylab="latitude", main="Point Xi DJF-97" )
  US(add=T)
  dev.off()