Let’s take a look at what the netcdf file looks like.

nc <- nc_open(nc_files[[1]])
nc 
## File /Users/dorh012/Documents/2020/hector-rcmip/output/rcmipII/netcdfs/rcmipII-1pctCO2-converted.nc (NC_FORMAT_NETCDF4):
## 
##      3 variables (excluding dimension variables):
##         float Surface Air Temperature Change[time,ensemble_member]   (Contiguous storage)  
##             _FillValue: NaN
##             climate_model: hector
##             region: World
##             _is_metadata: 0
##             unit: K
##             scenario: 1pctCO2
##             model: hector
##         float Surface Air Ocean Blended Temperature Change[time,ensemble_member]   (Contiguous storage)  
##             _FillValue: NaN
##             climate_model: hector
##             region: World
##             _is_metadata: 0
##             unit: K
##             scenario: 1pctCO2
##             model: hector
##         float Effective Radiative Forcing|Anthropogenic|Aerosols[time,ensemble_member]   (Contiguous storage)  
##             _FillValue: NaN
##             climate_model: hector
##             region: World
##             _is_metadata: 0
##             unit: W/m^2
##             scenario: 1pctCO2
##             model: hector
## 
##      2 dimensions:
##         time  Size:351
##             _FillValue: NaN
##         ensemble_member  Size:10000
##             units: ensemble_member
##             long_name: ensemble_member
## 
##     4 global attributes:
##         created_at: 2020-09-22 15:50:24
##         _scmdata_version: v1.2.0
##         scenario: 1pctCO2
##         model: hector
nc_close(nc)

Define a function that will extract and format the different data sets in prepration for plotting results.

lapply(nc_files, function(f){
  
  nc <- nc_open(f)
  
  ensemble <- as.character(0:9999)
  time     <- 1750:2100 # we know that this is going to be the time but we will only want to keep a subset of them
  yr_to_keep <- floor(seq(from=min(time), to=max(time), length.out = 50)) 
  
  vars <- c('Surface Air Temperature Change', 'Surface Air Ocean Blended Temperature Change', 
            'Effective Radiative Forcing|Anthropogenic|Aerosols')
  
  lapply(vars, function(v){
    
    var_output <- as.data.table(ncvar_get(nc, v))
    names(var_output) <- ensemble
    var_output[['time']] <- time
    var_output <- var_output[time %in% yr_to_keep, ]
    
    long <- melt(var_output, id.vars = 'time', variable.name = 'ensemble', value.name = 'value')
    long[['variable']] <- v
    
    return(long)
    
  }) %>%  
    dplyr::bind_rows()-> 
    data
  
  data$scenario <- ncatt_get(nc, 0)$scenario
  
  return(data)
  
  
}) %>% 
  dplyr::bind_rows() ->
  data
to_plot <- data[ , .("min" = min(value), "max" = max(value), "mean" = mean(value)), by = c("variable", "scenario", "time")]

Plots

The ribbon is the min/max hector results the dark black line is the Hector ensemble mean.

vars <- unique(to_plot$variable)
plot_list <- list()
for(v in vars){
  
  ggplot(data = to_plot[variable == v, ] ) + 
    geom_ribbon(aes(x = time, ymin = min, ymax = max),alpha = 0.5) +
    geom_line(aes(x = time, y=mean),size =2) +
    facet_wrap('scenario', scales= 'free') + 
    labs(title = paste0('RCMIP  II: ', v)) -> 
    plot
  
  plot_list[[v]] <- plot
  
}
plot_list$`Surface Air Temperature Change`

plot_list$`Surface Air Ocean Blended Temperature Change`

Let’s compare the two different types of temperature.

 ggplot(data = to_plot[variable %in% c('Surface Air Temperature Change',
                                       'Surface Air Ocean Blended Temperature Change' ), ] ) + 
    geom_ribbon(aes(x = time, ymin = min, ymax = max, fill = variable, color = variable),alpha = 0.5) +
    geom_line(aes(x = time, y=mean, color = variable)) +
    facet_wrap('scenario', scales = 'free') + 
    labs(title = paste0('RCMIP  II: ', v)) -> 
    plot
  plot

As we expected the mean surface temperature is cooler than the mean air temperature, and this difference is more pronounced in the future/higher \(CO_2\) scenarios.

plot_list$`Effective Radiative Forcing|Anthropogenic|Aerosols`

It makes snese that there is no uncertainty in the 1pctCO2, abrupt2xCO2 and the abrupt4xCO2 runs because they aren’t actually driven with aerosols.