# gfs ensemble forecast
# casestudy using GFS Ensembles
# Joseph wheatley 

library("ncdf")
library("sp")
library("rgdal")
library("maps")
library("maptools")
library("e1071")
library("diptest")
library("MASS")
library("lattice")

e <- c("01","02","03","04","05","06","07","08","09")
e <- c(e,paste(10:20))

day<-Sys.Date();
today <- paste(substr(day,1,4),substr(day,6,7),substr(day,9,10),sep="")
yesterday <- paste(substr(day,1,4),substr(day,6,7),as.numeric(substr(day,9,10))-1,sep="")

en.f <- function(i){
return(paste("ftp://ftpprd.ncep.noaa.gov/pub/data/nccf/com/gens/prod/gefs.",today,"/00/pgrb2a/gep",e[i],".t00z.pgrb2af384",sep=""))}

ens <- sapply(seq(1:20), en.f)

for(i in 1:20)
{

t<- "ensemble.grb"
tmp <- paste("ensemble",i-1,".nc",sep="")

loc=file.path(ens[i]);
download.file(loc,t,mode="wb"); # "wb" specifies binary file mode  
shell(paste("wgrib2 -s ",t," | grep \"TMP:2 m\" | wgrib2 -i ",t," -netcdf ",tmp,sep=""),intern=T)
}

t2m.en <- as.list(rep(0,20))
for (i in 1:20)
{

tmp <- paste("ensemble",i-1,".nc",sep="")
t2m <- open.ncdf(tmp)
t2m.en[[i]] <- get.var.ncdf(t2m,"TMP_2maboveground")
close.ncdf(t2m)
}

t2m <- open.ncdf(tmp)
x <- get.var.ncdf(t2m,"longitude")
y <- get.var.ncdf(t2m,"latitude")
close(t2m)



#array containing ensemble data
t2m.ens <- array(unlist(t2m.en),dim=c(360,181,20))-273

t2m_stats <- apply(t2m.ens,c(1,2), function(x) c(mean(x),median(x),sd(x),skewness(x,type=2),kurtosis(x,type=2),max(x),min(x),dip(x)))

Nstats <- dim(t2m_stats)[1]
# data frame longitude - latitude - ensemble member - sunlight - RH2m - T2m

######################
#  ensemble for some well known locations
##############################################

#change of coordinates degere E/W to deg E
shift <- function(c) { ifelse( c[1]<0,return(c+c(181,91)),return(c+c(1,91)))}
# cities long-lat
london <- c(0,51)
moscow <- c(38,56)
mumbai<- c(73,19)
beijing <- c(116,40)
rio <- c(-43,-23)
chicago <- c(-88,42)

interesting_places <- list(london,chicago,mumbai,rio,beijing,moscow)
placenames <- c("London","Chicago","Mumbai","Rio","Beijing","Moscow")
names(interesting_places) <- placenames
interesting_places <- lapply(interesting_places,shift)
#make a list 
place_ens <- as.data.frame(lapply( interesting_places, function(z) { t2m.ens[z[1],z[2],]}))

N_p <- length(placenames);

op=par(mfrow=c(3,2),oma=c(0,0,4,0))
for (i in 1:N_p){
hist(place_ens[[i]],col="firebrick",breaks=10,main=placenames[i],xlab="Temperature",ylab="Frequency",cex.main=1.5,font.lab=2,font.axis=2,cex.axis=1.4,cex.lab=1.3)
}
title(paste("2m Temperature 16 Day Ensemble\n", "GFS ", day, " 00UTC"),outer=T,cex.main=2.0,font.main=2)
par(op)


###########################
# convert to spatial
#########################
Nx <- length(x)
Ny <- length(y)
Nx2 <- Nx/2;
Nx2p <- Nx2+1;

t2m.sp <- array(rep(0,Nx*Ny*Nstats),dim=c(Nstats,Nx,Ny));
t2m.aux <- array(rep(0,Nx*Ny*Nstats),dim=c(Nstats,Nx,Ny));

t2m.aux <- t2m_stats[,,Ny:1];

t2m.sp[,1:Nx2,] <- t2m.aux[,Nx2p:Nx,];
t2m.sp[,Nx2p:Nx,] <- t2m.aux[,1:Nx2,];

x.c <- rep(0,Nx);
y.c <- rep(0,Ny);

# greenwich

x.c[1:Nx2] <- x[Nx2p:Nx]-360
x.c[Nx2p:Nx] <- x[1:Nx2]
y.c <- y

llCRS <- CRS("+proj=longlat +ellps=WGS84") 
gridPoints <- expand.grid(x.c,y.c)
sp <- SpatialPoints(gridPoints,proj4string = llCRS)

xmin <- min(x.c);
ymin <- min(y.c);

#  NOT satisfactory

c.offset <- c( -179.5,-89.5);
c.size <- c(1.0,0.994);
c.dim<-c(Nx,Ny);
grid <- GridTopology(c.offset,c.size,c.dim);

#grid <- coordinates(


dat <- data.frame("mean"=as.vector(t2m.sp[1,,]),"median"=as.vector(t2m.sp[2,,]),"sd"=as.vector(t2m.sp[3,,]),"skew"=as.vector(t2m.sp[4,,]),"kurtosis"=as.vector(t2m.sp[5,,]),"diptest"=as.vector(t2m.sp[8,,]));
sg <- SpatialGridDataFrame(grid, dat,proj4string = llCRS)

#use GDAL to produce a 

writeGDAL(sg, "ensembleIn.tif")
trans<- "C:/Progra~1/FWTools2.4.5/bin/gdalwarp -s_srs \"+proj=longlat +ellps=WGS84\" ensembleIn.tif ensembleStats.tif -t_srs \"+proj=cea +ellps=WGS84\""
shell(trans)

ll <- readGDAL("ensembleStats.tif")
proj4string(ll) <- CRS("+proj=cea +ellps=WGS84")



wrld <- map("world",xlim=c(-179,179),ylim = c(-89,89),plot=F)
wrld_p <- pruneMap(wrld,xlim =c(-179,179))
wrld_sp <- map2SpatialLines(wrld_p,proj4string = llCRS)
land <- list("sp.lines", wrld_sp, first=FALSE, col="darkgreen")

wrld_cylin <- spTransform(wrld_sp,CRS("+proj=cea +ellps=WGS84"))
land_cylin <- list("sp.lines", wrld_cylin, first=FALSE, col="darkgreen")
land_cylin1 <- list("sp.lines", wrld_cylin, first=FALSE, col="darkgreen",alpha=1.0)

gr <- gridlines(wrld_sp,easts=seq(-180,180,30),norths=seq(-90,90,30))
gr_cylin<- list("sp.lines",spTransform(gr,CRS("+proj=cea +ellps=WGS84")),first=F,col="gray1",lty=2)

x.ll <- seq(from=ll@grid@cellcentre.offset[1], by = ll@grid@cellsize[1], length.out = ll@grid@cells.dim[1])
y.ll <- seq(from=ll@grid@cellcentre.offset[2], by = ll@grid@cellsize[2], length.out = ll@grid@cells.dim[2])

contourdat <- matrix(ll@data$band2,ll@grid@cells.dim[1],ll@grid@cells.dim[2]);
contourdat <- contourdat[,ll@grid@cells.dim[2]:1]

cc <- contourLines(x.ll,y.ll,contourdat,levels=seq(-10,30,by=5))
cc <- ContourLines2SLDF(cc, proj4string=CRS("+proj=cea +ellps=WGS84"))
cc_labs <- toString(seq(from=240,to=310,by=5))

ll@bbox[2,] <- wrld_cylin@bbox[2,]
cc@bbox[2,] <- wrld_cylin@bbox[2,]

cc_cylin <- list("sp.lines",cc,first=F,col="gold",lwd=1,alpha=0.2)
sd.palette <- colorRampPalette(c("black","orange","firebrick"), space = "rgb")
#skew.palette <- colorRampPalette(c("violet","black","red"), space = "rgb")
skew.palette <- brewer.pal(9,"PuOr")

s1 <- spplot(ll,c("band3"),sp.layout=list(land_cylin,gr_cylin,cc_cylin),main=paste("ensemble standard deviation"),col.regions=sd.palette,cuts=200)
s2 <- spplot(ll,c("band4"),sp.layout=list(land_cylin1,gr_cylin,cc_cylin),main=paste("ensemble skew"),col.regions=skew.palette,cuts=7)

print(s1, c(0,0.5,1,1),more=T)
print(s2,c(0,0,1,0.5))