NCEP Global Forecast System

Just about everyone is familiar with weather maps. There are many situations where it is useful to combine the underlying numerical weather data with other types of information. Accessing  the weather data is a necessary first step.

The output from the U.S. National Centers for Environmental Prediction (NCEP) Global  Forecast System (GFS) is freely available. The surface resolution of the model is ≈ 0.3º× 0.3°. The model runs every 6 hours, producing forecasts at 3-hourly intervals extending out to 16 days. As an example of output from GFS, the map (below) shows the predicted average  temperature at 2 metres over the entire globe for the next 24 hr (date of this post). The map shows predicted cold conditions in Europe, and the continuing heatwave in Australia.

t2m

How the map was made

GFS forecasts are in a format called GRIB2. According to Wikipedia, “GRIB (GRIdded Binary) is a mathematically concise data format commonly used in meteorology to store historical and forecast weather data.” GRIB files contain physical fields such as temperature, humidity etc defined on a spatial grid, as well as boundary conditions such as vegetation type and elevation. The data might be assimilated from observations, or output from a forecast model.

The first step is to translate the GRIB into a raster format such as netcdf which can be read in R. For example, the GRIB2 file gfs.2009121700/gfs.t00z.sfluxgrbf03.grib2 contains the 3-hr forecast surface data on 17 Dec 2009  produced at 00 UTC (midnight universal time). An inventory of the data contained in this file can be seen here. Download this forecast as temp.grb

loc=file.path("ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.2009121700/gfs.t00z.sfluxgrbf03.grib2")
download.file(loc,"temp.grb",mode="wb")

To read temp.grb a utility called wgrib2 needs to be installed on your system. Then data such as land fraction can extracted into a netcdf file LAND.nc using the R shell command

shell("wgrib2 -s temp03.grb | grep :LAND: | wgrib2 -i temp00.grb -netcdf LAND.nc",intern=T)

The ncdf package can now be used to read the contents of LAND.nc.

library(ncdf)
landFrac <-open.ncdf("LAND.nc")
land <- get.var.ncdf(landFrac,"LAND_surface")
x <- get.var.ncdf(landFrac,"longitude")
y <- get.var.ncdf(landFrac,"latitude")

The 1152×576 matrix land takes values 1 for land and 0 for water (sea-ice is 1). x and y are the longitude and latitude of the non-uniform GFS grid.

2m temperature data can be read in the same way. The average of the first 8  forecasts was called t2m.mean and plotted using image.plot() from the fields package:

library(fields)
rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")#colors
image.plot(x,y,t2m.mean,col=rgb.palette(200),axes=F,main=as.expression(paste("GFS 24hr Average 2M Temperature",day,"00 UTC",sep="")),axes=F,legend.lab="o C")
contour(x,y,land,add=TRUE,lwd=1,levels=0.99,drawlabels=FALSE,col="grey30") #add land outline

December 16, 2009  Tags: , ,   Posted in: Climate, Uncategorized  4 Comments

How big is the uncertainty in the global temperature trend?

The global temperature record shows strong serial correlation. This boosts the uncertainty in the slope of the global warming trend. Here the uncertainty in the trend is derived using one of the satellite datasets.

warming in the satellite era

Data from satellites have been used to construct a global temperature record beginning in December 1978. This covers the recent warming period. The University of Alabama at Huntsville (UAH) 2LT product gives average temperature of the lower atmosphere derived from microwave radiances. Planck’s law says that microwave radiances are related to temperature.

The plot below shows UAH monthly time-series and a least-squares regression fit. The slope of the regression line is 0.013ºC/year (equivalent to 0.13°C/decade or 1.3°C/century). Also shown is the Gistemp global temperature time-series over the same time interval. Gistemp is based on surface weather stations. The trend in the Gistemp series is higher, 0.18°C/decade.

anomaly

uncertainty in trend

A linear regression fit to global temperatures is a statistical model in which strong residual variations are superposed on a linear warming trend. The residual variations reflect “noise” in the climate system. ENSO, random events such as volcanos, and many other modes of variability contribute to the “noise”. The “uncertainty” in the global warming trend can be estimated in the context of this statistical model.

From this point of view, the UAH monthly data are one realisation of a random sample of length N. In reality we only have one sample, the UAH data from the real world. However, if repeated sampling and slope measurement were possible, it would yield a range of outcomes due to noise. The width of the resulting slope probability distribution reflects the uncertainty in the global warming trend. It is evident that there is strong serial correlation in the residuals. This enhances the uncertainty, because it reduces the effective number of independent degrees of freedom in the data. Fewer degrees of freedom means more uncertainty.

The UAH residuals are fit to an auto-regressive model to account for serial auto-correlation. The Akaike information criterion to choose an optimal order n. For UAH it  turns out that n=2, with coefficients 0.58  and 0.24. To find the slope probability distribution 10,000 time-series were used in a Monte Carlo simulation. Each simulated time-series assumes a UAH trend-line of 0.13°C/decade and an AR(2) noise term parameterized using the AR(2) model fit to the UAH residuals. Linear regression fit on each of these simulated time-series yields the slope probablity distribution.

slope

A typical result is shown above along with a fit to the normal distribution. By construction the mean is 0.13°C/decade as expected. The standard error is 0.036°C/decade. This is 3.4 times greater than the standard error with independent and identically distributed (iid) residuals.

The 2σ confidence interval for UAH global temperature trend is 0.05 to 0.2°C/decade. Unfortunately, this is quite a wide window. Would another 30 years of data help? Assuming that the properties of the UAH time-series do not change, and extending the above simulation for another 30 years, the standard error is reduced from 0.035°C/decade to 0.026°C/decade. Better, but still a surprising degree of uncertainty.

R code

The above analysis was carried out using R functions ar() and arima.sim()

fitLin <-lm(GLOBAL~dates) # linear regression fit to GLOBAL= UAH global temperature, dates = UAH dates
ar.fit <- ar(residuals(fitLin)) # AR model for residuals
N=10000; # number of simulated UAH time-series
fits <- cbind(rep(0,N),rep(0,N))
trend <- coefficients(fitLin)[1] + dates*coefficients(fitLin)[2] # UAH trend line
#Monte Carlo
for( i in 1:N){
warming <- trend+arima.sim(n = NData, list(ar = c(ar.fit$ar[1], ar.fit$ar[2])), sd=sqrt(ar.fit$var.pred)) #generate sample
fit.warm <- lm(warming~dates) # linear fit
fits[i,1] <- coefficients(fit.warm)[1];
fits[i,2] <- coefficients(fit.warm)[2]; # slope of simulated data
}
library(MASS) # MASS library must be installed
gf <- fitdistr(fits[,2],"normal") # fit a normal distribution to the slopes
mn <-gf$estimate[1] #mean
sig <- gf$estimate[2] #standard error

December 8, 2009   Posted in: Climate, Uncategorized  No Comments