{\rtf1\ansi\ansicpg1252\cocoartf1038\cocoasubrtf360
{\fonttbl\f0\fmodern\fcharset0 Courier;}
{\colortbl;\red255\green255\blue255;}
\paperw11900\paperh16840\margl1440\margr1440\vieww16540\viewh14300\viewkind0
\deftab720
\pard\pardeftab720\ql\qnatural

\f0\fs24 \cf0 # 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  \
system(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,paste("TMP__",i,"_2maboveground",sep=""))\
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))\
\
\
}