Verifying Ensembles

Read...

Once you have an ensemble, you need to evaluate its skill. Two methods of verifying ensembles will be discussed in this lesson: the rank histogram and the continuous ranked probability score. Both methods evaluate the spread of the ensemble and can be used to assess the skill of the forecast. Read on to learn more.

Rank Histogram Method

A rank histogram is a diagnostic tool to evaluate the spread of an ensemble. For each forecast, bins are selected by ranking the ensemble member forecast (lowest to highest). For N ensemble members, there will be N+1 bins - the lowest and highest bins are open-ended. (source: The Rank Histogram) Then you create a frequency histogram by plotting the number of occurrences in each bin.

A rank histogram is a great plot to eyeball if the ensemble is over or under-dispersed, or to see if the ensemble is biased. In R, you can use the package ‘SpecsVerification’ to create a rank histogram.

Below are several different examples of rank histograms. Match them to the most appropriate description. 

Continuous Ranked Probability Score (CRPS)

A Continuous Ranked Probability Score (CRPS) is the PDF equivalent of the MSE. The CRPS rewards forecasts based on the amount of probability at an observed value. The higher the density, the lower the score. CRPS indicates how close the observation is to the weighted average of the PDF. A lower CRPS indicates better model performance (the PDF matches the distribution of forecast errors) – thus, CRPS lets you compare different ensemble models to pick the best. In R, you can use the ‘SpecsVerification’ package to create this plot.

For more information, I suggest you check out this link.

Example

Let’s create these two graphs in R. We will start with the rank histogram. Begin by downloading this GEFS forecast output for NYC. The file contains 24-hour forecast temperatures valid at 00Z from January 1st through November 30th, 2018. Start by loading in the file.

# load in temperature forecast from GFS for NYC
library(RNetCDF)
fid <- open.nc('tmp_2m_latlon_all_20180101_20181130_jzr8iFlxW9.nc')
mydata <- read.nc(fid)

Now, extract out the forecast valid time.

validTime <- mydata$intValidTime

We want to match up the observations at the time for which the forecast was valid. First, download this data set of hourly temperature observations. Load in the data set and convert the time variable to match that of the valid time forecast variable (YearMonthDayHour) by using the code below:

# load in daily temperature observations
mydata2 <- read.csv("NCYTemps2018.csv")
obsTime <- as.POSIXct(strptime(mydata2$Date,"%m/%d/%Y %H:%M"))
tempTime <- format(obsTime,"%Y%m%d%H")

Now, loop through each forecast and pull out the corresponding observation that matches the time the forecast is valid.

# for each forecast valid time find the corresponding observation time
obsTemp <- array(NA,length(validTime))
for(itime in 1:length(validTime)){
  goodIndex <- which(as.character(validTime[itime]) == tempTime)
  obsTemp[itime] <- mydata2$AirTemp[goodIndex]
}

Next, convert the forecast temperature to Fahrenheit (from Kelvin) and transpose the matrix, so the format is (n,k) where n is the forecast and k is the ensemble member. This is the setup needed to use the rank histogram function in R.

# convert kelvin to F
forecastTemp <- (9/5)*(mydata$Temperature_height_above_ground-273.15)+32
# transpose forecast temp (n,k) where k is the number of members
forecastTemp <- t(forecastTemp)

Now, we can create the rank histogram. We will first estimate the histogram using the ‘Rankhist’ function, and then plot the histogram using the ‘PlotRankhist’ function. Run the code below:

# load in library
library('SpecsVerification')
# create rank histogram
rh <- Rankhist(forecastTemp,obsTemp,handle.na="use.complete")
PlotRankhist(rh)

You should get the following figure:

Rank histogram for GEFS forecasted temperature in NYC for January 1st through November 30th, 2018 valid at 00Z

Credit: J. Roman

This is an example of an ensemble that has too small of a spread.

Now, let’s make a CRPS plot. Begin by downloading this forecast output of temperature for NYC. Again, the forecast output is from January 1st, 2018 – November 30^th, 2018 but, this time, we have the output for all forecast times from 0 hours out to 192 hours (8 days). We will use the same observational data set as before. Begin by loading in the temperature forecast.

# load in temperature forecast from GFS for NYC
library(RNetCDF)
fid <- open.nc('tmp_2m_latlon_all_20180101_20181130_jzr8YY8Qhv.nc')
mydata <- read.nc(fid)

Now, extract the time and select out the forecast hours you are interested in testing. I’m going to look at the 1-day forecast, 2-day, 3-day, and so on to the 8-day forecast. Use the code below to set this up and initialize the CRPS variable.

time <- as.Date(mydata$time/24,origin="1800-01-01 00:00:00")
fhour <- seq(24,192,24)
tempTime <- format(obsTime,"%Y%m%d%H")

# initialize variable
CRPS <- array(NA,length(fhour))

Now, we are ready to loop through each forecast hour and estimate the mean CRPS. Use the code below:

# loop through each fhour and estimate the mean CRPS
for(ihour in 1:length(fhour)){
  fhourIndex <- which(mydata$fhour == fhour[ihour])
  # extract out valid times
  validTime <- mydata$intValidTime[fhourIndex,]
  # match up observations
  obsTemp <- array(NA,length(time))
  for(itime in 1:length(validTime)){
    goodIndex <- which(as.character(validTime[itime]) == tempTime)
    if(length(goodIndex) == 0){
      obsTemp[itime] <- NA
    }else{
      obsTemp[itime] <- mydata2$AirTemp[goodIndex]
    }
    
  }
  
  # convert kelvin to F
  forecastTemp <- (9/5)*(mydata$Temperature_height_above_ground[fhourIndex,,]-273.15)+32
  # transpose forecast temp (n,k) where k is the number of members
  forecastTemp <- t(forecastTemp)
  # compute the CRPS
  tempCRPS <- EnsCrps(forecastTemp,obsTemp)
  CRPS[ihour] <- mean(tempCRPS,na.rm=TRUE)
}

The code starts by pulling out the valid times for the particular forecast hour of interest. Then we match up the observations for the valid forecast time. Next, we convert the forecast temperature to Fahrenheit and set up the matrix for the forecast values to be in the format (n,k) where n is the forecast and k is the ensemble member. Finally, you compute the CRPS using the function ‘EnsCrps’ and save the mean value to our CRPS variable.

Now, plot the mean CRPS for each forecast hour using the code below:

# plot
plot.new()
grid()
par(new=TRUE)
plot(fhour/24,CRPS,xlim=c(1,8),ylim=c(3,4.5),type="b",pch=20,lty=1,lwd=2,
     xlab="Forecast Day",ylab="CRPS (degrees)",
     main="GEFS Temperature Forecast NYC")

You should get the following figure:

CRPS plot for GEFS forecasted temperature in NYC from January 1st through November 30th, 2018 for various forecast days

Credit: J. Roman

As a reminder, a lower CRPS value is better. Thus, the 1-day forecast is best (as expected) and the 8-day forecast is the worst (again as expected). By adding in more ensemble models (such as the ECMWF) we could assess which model is best, based on which one has the lowest CRPS value, for each forecast hour.