R and GFS.

Playing with NOAA's GFS data with R.

Credit:

• http://braz.blogspot.com/2009/12/r-for-global-climate-forecasting.html 
• http://joewheatley.net/ncep-global-forecast-system/

Here are some sample images.

Pre-requisites

• Fetched data from NOAA (in particular this dataset).
• Installed wgrib2.
• Extracted the fields I was interested in from the grib2 file into netcdf

./grib2/wgrib2/wgrib2  -s  gfs.t12z.sfluxgrbf03.grib2  | grep :LAND | ./grib2/wgrib2/wgrib2 -i gfs.t12z.sfluxgrbf03.grib2 -netcdf land.netcdf
./grib2/wgrib2/wgrib2  -s  gfs.t12z.sfluxgrbf03.grib2  | grep TMP:2 | ./grib2/wgrib2/wgrib2 -i gfs.t12z.sfluxgrbf03.grib2 -netcdf temp2.netcdf
./grib2/wgrib2/wgrib2  -s  gfs.t12z.sfluxgrbf03.grib2  | grep VEG: | ./grib2/wgrib2/wgrib2 -i gfs.t12z.sfluxgrbf03.grib2 -netcdf vegetation.netcdf
./grib2/wgrib2/wgrib2  -s  gfs.t12z.sfluxgrbf03.grib2  | grep TCDC: | ./grib2/wgrib2/wgrib2 -i gfs.t12z.sfluxgrbf03.grib2 -netcdf cloud.netcdf

• You can also view the available fields using:

./grib2/wgrib2/wgrib2  -s  gfs.t12z.sfluxgrbf03.grib2  | less

• Or by visiting the description.

Having fun

Started R and installed netcdf and fields:

install.packages('ncdf')
install.packages('fields')

Then used the following gist in R.

	# load datasets
	library(fields)
	library(ncdf)
	# Load land data
	landFrac <-open.ncdf("land.netcdf")
	land <- get.var.ncdf(landFrac,"LAND_surface")
	x <- get.var.ncdf(landFrac,"longitude")
	y <- get.var.ncdf(landFrac,"latitude")
	# Load temp data
	tempFrac <-open.ncdf("temp2.netcdf")
	# Temperature 2 metres above ground
	temp_k <- get.var.ncdf(tempFrac, "TMP_2maboveground")
	temp_c <- temp_k - 273.15 # Kelvins to Centigrade
	# Palette for temp data
	rgb.palette.temp <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")
	# View temp plot
	image.plot(x,y,temp_c,col=rgb.palette.temp(200),axes=F,main=as.expression(paste("GFS 24hr Average 2M Temperature ",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="o C")
	# save plot to file
	png("/tmp/temp.png", 1600, 800)
	image.plot(x,y,temp_c,col=rgb.palette.temp(200),axes=F,main=as.expression(paste("GFS 24hr Average 2M Temperature ",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="o C")
	dev.off()
	# Vegetation
	vegFrac <-open.ncdf("vegetation.netcdf")
	vegetation <- get.var.ncdf(vegFrac, "VEG_surface")
	# Truncate unrealistic values
	vegetation[vegetation > 100000] <- 0
	# Palette
	rgb.palette.veg <- colorRampPalette(c("seashell","firebrick","orange","seagreen"), space = "rgb")
	# Plot
	image.plot(x,y,vegetation,col=rgb.palette.veg(100),axes=F,main=as.expression(paste("Percentage of land under vegetation cover",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="coverage %")
	# Or save png
	png("/tmp/vegetation.png", 1600, 800)
	image.plot(x,y,vegetation,col=rgb.palette.veg(100),axes=F,main=as.expression(paste("Percentage of land under vegetation cover",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="coverage %")
	dev.off()
	# Cloudiness
	cloudFrac <-open.ncdf("cloud.netcdf")
	lowclouds <- get.var.ncdf(cloudFrac, "TCDC_lowcloudlayer")
	rgb.palette.clouds <- colorRampPalette(c("royalblue1", "lightskyblue", "lightsteelblue1", "papayawhip","snow"), space = "rgb")
	# Preview
	image.plot(x,y, lowclouds, col=rgb.palette.veg(100),axes=F,main=as.expression(paste("Percentage of total cloud cover (low clouds) ",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="coverage %")
	image.plot(x,y, lowclouds, col=rgb.palette.clouds(100),axes=F,main=as.expression(paste("Percentage of total cloud cover (low clouds) ",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="coverage %")
	contour(x,y,land,add=TRUE,lwd=1,levels=0.99,drawlabels=FALSE,col="grey30") #add land outlineland2 <- get.var.ncdf(landFrac2, "LAND_surface")
	# Save png
	png("/tmp/lowclouds.png", 1600, 800)
	image.plot(x,y, lowclouds, col=rgb.palette.clouds(100),axes=F,main=as.expression(paste("Percentage of total cloud cover (low clouds) ",format(Sys.time(), "%Y-%m-%d") ," 03 UTC",sep="")), legend.lab="coverage %")
	contour(x,y,land,add=TRUE,lwd=1,levels=0.99,drawlabels=FALSE,col="grey30") #add land outlineland2 <- get.var.ncdf(landFrac2, "LAND_surface")
	dev.off()

view raw noaa_gfs.R hosted with ❤ by GitHub