#usage:  R BATCH --vanilla reanbin.R
#reanbin.R will read in a binary file, in the size and form of the surface (gaussian)
#reanalysis grid.  If the binary file came from GrADS, via the "set gxout fwrite",
#note the GrADS adds a redundant "column" in the extreme east, so the size is 193*94.
#Reanalyis data directly from wgrib is 192*94.
#GrADS flips the data north-south, from wgrib.
#Note the companion source code "gaussgrid.R" is required.
#The contour plotting of R is very inefficient in time and storage;
#Here the huge .eps file is converted to a .png and the .eps is deleted.
#Obviously, this script could be extended to do some more statistical processing.
#Only simple plotting is done here, for which GrADS would have been quite capable of.
#wgrib example:
# wgrib /REANALYSIS/1997/data/at00z12z/tmp.2m -d 400 -V -o fromwgrib.dat -nh
library("fields") #needed for map overlay
source("gaussgrid.R")
#
#use the following two lines for data from GrADS
nlon=193
binfilename='climo7.dat'
#
#use the following two lines for data from wgrib
#nlon=192
#binfilename='fromwgrib.dat'
#
nlat=94
nm=nlon*nlat
s=readBin(binfilename,numeric(),n=nm,size=4)
#matrix storage convention is rotated 90 degrees from map: 
gdat=matrix(s,nrow=nlon,ncol=nlat)
#inspect the data:
str(gdat)
max(gdat)
min(gdat)
gdat[nlon,45]
gdat[1,45]
thedat=gdat
if (nlon==192) thedat=gdat[1:nlon,94:1] #flip over if direct from wgrib
thelats<-gausslats[94:1] #must flip over the latitudes for both types of files
thelons=gausslons[1:nlon]
tlevs=seq(170,320,10)
ncolors=length(tlevs)-1
somecolors=rainbow(ncolors,start=0.,end=0.7)
thecolors=somecolors[ncolors:1] #reverse colors, so high values are red
postscript(file="batch.eps",onefile=FALSE,horizontal=FALSE,width=6.0,height=4.0)
filled.contour(thelons,thelats,thedat,xlab='lon',ylab='lat',
  col=thecolors,levels=tlevs,
  plot.title={title(binfilename)},
  plot.axes={
    world(ylim=c(-90,90), xlim=c(0,360), add=TRUE, 
      xlab = "", ylab = "",  xaxt = "n", yaxt = "n", 
      eps=.1,shift=TRUE); axis(1); axis(2)
  }
)
dev.off()
system("convert +antialias -density 144x144 batch.eps batch.png")
system("rm batch.eps")