#############################################################################################
# jet.R produces graphics of today's jet stream and two day forecast.
# Based on NCEP GFS Analysis at 00 UTC
# FAST DOWNLOAD using get_gfs.pl http://www.cpc.noaa.gov/products/wesley/get_gfs.html (w. ebisuzaki)
# For use of sp package, see Bivand, Pedesma & Gomez-Rubio Applied Spatial Data Analysis with R
# Joe Wheatley Biospherica
############################################################################################

library(ncdf)
library(sp)
library(maps)
library(maptools)
library(rgdal)


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="")

llCRS <- CRS("+proj=longlat +ellps=WGS84") #coordinate reference system

# northern hemisphere country map & grid

hemi_north <- map("world",xlim=c(-179,179),ylim = c(0,89),plot=F)
hemi_north_p <- pruneMap(hemi_north,xlim =c(-179,179))
hemi_north_sp <- map2SpatialLines(hemi_north_p,proj4string = llCRS)

stereo_north <- CRS("+proj=ups +north")  # universal polar stereographic
polar_north <- spTransform(hemi_north_sp,stereo_north)

# add grid lines
gr_sp_north <- gridlines(hemi_north_sp,easts=seq(-180,180,20),norths=seq(0,80,20))
gr_polar_north <- spTransform(gr_sp_north,stereo_north)

# southern hemisphere country map & grid

hemi_south <- map("world",xlim=c(-179,179),ylim = c(-89,0),plot=F)
hemi_south_p <- pruneMap(hemi_south,xlim =c(-179,179))
hemi_south_sp <- map2SpatialLines(hemi_south_p,proj4string = llCRS)

stereo_south <- CRS("+proj=ups +south")  # universal polar stereographic
polar_south <- spTransform(hemi_south_sp,stereo_south)

gr_sp_south <- gridlines(hemi_south_sp,easts=seq(-180,180,20),norths=seq(-80,0,20))
gr_polar_south <- spTransform(gr_sp_south,stereo_south)

# download 300mb wind velocity components for "00" today 00UTC "24" tomorrow 00UTC etc

#forecast hours
fc = c("00","24","48"); 
#perl script get_gfs.pl is located in current directory 
getgfs<- paste("C:/strawberry/perl/bin/perl get_gfs.pl data ", today,"00 0 48 24 UGRD:VGRD 300_mb  jet_gfs",sep="")  
shell(getgfs,intern=T)

# netcdf wind component files
for(i in 1:length(fc))
{
t <- paste("jet_gfs/gfs.t00z.pgrb2f",fc[i],sep=""); #grib2 files

shell(paste("wgrib2 -s ",t," | grep UGRD:300 | wgrib2 -i ",t," -netcdf Ucomp",fc[i],".nc",sep=""),intern=T) #decode grib2 & extract zonal wind speed Uvel
shell(paste("wgrib2 -s ",t," | grep VGRD:300 | wgrib2 -i ",t," -netcdf Vcomp",fc[i],".nc",sep=""),intern=T) #decode grib2 & extract meridonal wind spped Vvel
#shell(paste("wgrib2 -s ",t," | grep TMP:300 | wgrib2 -i ",t," -netcdf T300.nc",sep=""),intern=T) #decode grib2 & extract meridonal wind spped Vvel
}

jet.palette <- colorRampPalette(c("lavender","steelblue","purple"), space = "rgb")

arctic <- as.list(rep(0,length(fc)))
antarctic <-as.list(rep(0,length(fc)))

for( i in 1:length(fc))
{
Uc <- open.ncdf(paste("Ucomp",fc[i],".nc",sep=""));
Vc <- open.ncdf(paste("Vcomp",fc[i],".nc",sep=""));

Uvel <- get.var.ncdf(Uc,"UGRD_300mb") #2m temperature maplist 
Vvel <- get.var.ncdf(Vc,"VGRD_300mb")

x_jet <- get.var.ncdf(Uc,"longitude")
y_jet <- get.var.ncdf(Uc,"latitude")

close.ncdf(Uc);
close.ncdf(Vc);

Nx_jet <- length(x_jet);
Ny_jet <- length(y_jet);

i_r <- seq(from=1, to=Nx_jet,by=10);
j_r <- seq(from=1,to=Ny_jet,by=10);

xy.g <- expand.grid(x_jet[i_r],y_jet[j_r]);
jet <- sqrt(Uvel^2+Vvel^2); # absolute velocity

jet_decimated <- jet[i_r,j_r];

# creat sp jetstream object

jet <- sqrt(Uvel^2+Vvel^2);
jet <- jet[,Ny_jet:1]
Nx2 <- Nx_jet/2;
Nx2p <- Nx2+1;
jet1 <- matrix(rep(0,Nx_jet*Ny_jet),Nx_jet,Ny_jet);
jet1[1:Nx2,] <- jet[Nx2p:Nx_jet,];
jet1[Nx2p:Nx_jet,] <- jet[1:Nx2,];

x.c <- rep(0,Nx_jet);
y.c <- rep(0,Ny_jet);

# greenwich

x.c[1:Nx2] <- x_jet[Nx2p:Nx_jet]-360
x.c[Nx2p:Nx_jet] <- x_jet[1:Nx2]
y.c <- y_jet

gridPoints <- expand.grid(x.c,y.c)
sp <- SpatialPoints(gridPoints,proj4string = llCRS)

xmin_jet <- min(x.c);
ymin_jet <- min(y.c);

#############  NOT satisfactory

c.offset <- c( -179.75,-89.75);
c.size <- c(0.5,0.4985);
c.dim<-c(Nx_jet,Ny_jet);
grid <- GridTopology(c.offset,c.size,c.dim);

dat <- data.frame("band1"=as.vector(jet1));
sg <- SpatialGridDataFrame(grid, dat,proj4string = llCRS)


#south 

sg_south <- spTransform(sg, CRS("+proj=ups +south"))
sg_south@bbox <- polar_south@bbox

#north

sg_north <- spTransform(sg, CRS("+proj=ups +north +units=m"))
sg_north@bbox <- polar_north@bbox


# build a SpatialGridDataFrame from sg_polar using akima linear interpolation

library(akima)

# north

N_res <- 200 #resolution of new grid

reg <- bbox(sg_north)
x_in <- sg_north@coords[,1]
y_in <- sg_north@coords[,2]
dat_in <- sg_north@data$band1

x_in_dec <- x_in[seq(1,length(x_in),by=4)]
y_in_dec <- y_in[seq(1,length(y_in),by=4)]
dat_in_dec <- dat_in[seq(1,length(dat_in),by=4)];


#fit <- Tps(as.matrix(pts),dats,df=500) # thin plate spline
#out <- predict.surface(fit,nx=400,ny=400)

x0 <- seq(from=reg[1,1],to=reg[1,2],length.out=N_res)
y0 <- seq(from=reg[2,1],to=reg[2,2],length.out=N_res)

int <- interp(x_in_dec,y_in_dec,dat_in_dec,x0,y0)
#image(x0,y0,int$z)


newgrid_north <- GridTopology(c(reg[1,1],reg[2,1]),c((reg[1,2]-reg[1,1])/N_res,(reg[2,2]-reg[2,1])/N_res),c(N_res,N_res))
newdata_north <- data.frame("band1"= as.vector(int$z[,N_res:1]))  # flip y-axis required for some reason
new_sp_north <- SpatialGridDataFrame(newgrid_north, newdata_north,proj4string = CRS("+proj=ups +north +units=m") )

land_north <- list("sp.lines", polar_north, first=FALSE,col="darkgreen")
gridlines_north <- list("sp.lines", gr_polar_north,first=FALSE,col="gray1",lty=2)

# south


reg <- bbox(polar_south)
x_in <- sg_south@coords[,1]
y_in <- sg_south@coords[,2]
dat_in <- sg_south@data$band1

x_in_dec <- x_in[seq(1,length(x_in),by=4)]
y_in_dec <- y_in[seq(1,length(y_in),by=4)]
dat_in_dec <- dat_in[seq(1,length(dat_in),by=4)];


#fit <- Tps(as.matrix(pts),dats,df=500) # thin plate spline
#out <- predict.surface(fit,nx=400,ny=400)

x0 <- seq(from=reg[1,1],to=reg[1,2],length.out=N_res)
y0 <- seq(from=reg[2,1],to=reg[2,2],length.out=N_res)

int <- interp(x_in_dec,y_in_dec,dat_in_dec,x0,y0)
#image(x0,y0,int$z)

newgrid_south <- GridTopology(c(reg[1,1],reg[2,1]),c((reg[1,2]-reg[1,1])/N_res,(reg[2,2]-reg[2,1])/N_res),c(N_res,N_res))
newdata_south <- data.frame("band1"= as.vector(int$z[,N_res:1]))  # flip y-axis required for some reason
new_sp_south <- SpatialGridDataFrame(newgrid_south, newdata_south,proj4string = CRS("+proj=ups +south +units=m") )

land_south <- list("sp.lines", polar_south, first=FALSE, col="darkgreen")
gridlines_south <- list("sp.lines", gr_polar_south,first=FALSE,col="gray1",lty=2)

#plot 300mb winds for both hemispheres

arctic[[i]]<-spplot(new_sp_north,col.regions=jet.palette, cuts=100,at=seq(0,120,by=1),sp.layout = list(land_north,gridlines_north),main=paste(day+i-1," 00UTC",sep=""))
antarctic[[i]]<-spplot(new_sp_south,col.regions=jet.palette, cuts=100,at=seq(0,120,by=1),sp.layout = list(land_south,gridlines_south),main=paste(day+i-1," 00UTC",sep=""))

}
#make table of plots
Nt <- length(fc);
Ntm <- Nt-1;
for(i in 1:Ntm)
{
print(arctic[[i]],c(0,1-i/Nt,0.6,1-(i-1)/Nt),more=T)
print(antarctic[[i]],c(0.4,1-i/Nt,1,1-(i-1)/Nt),more=T)
}

print(arctic[[Nt]],c(0,0,0.6,1-(Nt-1)/Nt),more=T)
print(antarctic[[Nt]],c(0.4,0,1,1-(Nt-1)/Nt))