Plotting TraCE-21k climate map

Information

In this document you will find information on how to extract data from a TraCE-21ka file and plot a series of climate maps for a specific time period (5000-7000 cal BP). This guideline will use the example of surface temperature, however this can easily be changed for a different climate variable and or model file.

You can download the netCDF file here. You will need to make a free online account to download. The file names are relatively long, but the one you are looking for is: “trace.01-36.22000BP.cam2.TS.22000BP_decavg_400BCE.nc” for temperature… where it says ‘TS’ in the center of this file name is where it refers to temperature, so be sure you have this one!

This code is built in Quartro to visually show both the code and output, however to use this yourself simply copy the code snippets and paste in your own RStudio script.

If you happen to use parts of this code, please ensure to credit it! Any questions please email at laura.boyall.2016@live.rhul.ac.uk

Happy plotting!

Loading the required packages

R requires a series of packages to be loaded to help with streamlining the different functions. The following packages are the ones which I find useful, but if you are more familiar with alternatives, then please do use.

You only have to install packages once, however each time you begin a new R session, you have to retrieve them from the library using the library(package name) function.

Remember to then set your working directory as normal to where your netCDF files are stored.

if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("ncdf4", quietly = TRUE)) install.packages("ncdf4")
if (!requireNamespace("raster", quietly = TRUE)) install.packages("raster")
if (!requireNamespace("dplyr", quietly = TRUE)) install.packages("dplyr")
if (!requireNamespace("maps", quietly = TRUE)) install.packages("maps")
if (!requireNamespace("maptools", quietly = TRUE)) install.packages("maps")
if (!requireNamespace("sf", quietly = TRUE)) install.packages("maps")
if (!requireNamespace("rgdal", quietly = TRUE)) install.packages("maps")

library(ggplot2) # used for plotting
library(ncdf4) # used for reading and extracting data from NetCDF model file
library(raster) # helpful package for plotting maps
library(dplyr) # helpful tool for filtering data
library(maps) # to create the map outlines
library(sf) # essential for spatial data handing
library(maptools) # converting data to polygons
library(rgdal) # used to help with map projections

Reading in the data

The data can be downloaded from here. Whilst monthly data is available, the file sizes here can be rather large, so decadal means should be sufficient. Note that there are atmospheric, land, oceanic and ice simulations… if you cannot find the variable of interest in one, check the other! (Temperature (TS), precipitation (PRECL), sea level pressure (SLP) etc is stored in atmospheric).

If you are using this on your own machine, change the source directory to where your model file is stored.

Temperature_data <- nc_open("trace.01-36.22000BP.cam2.TS.22000BP_decavg_400BCE.nc")

Now that we have uploaded our model file into R, we can extract the variable of interest. Each model file is going to be completely different. Some will only have one variable to plot e.g. temperature, however this still may have a number of different names. Some examples I have come across is TS, TAS, tas, T, temp, TC. A good way of checking what your variable is stored as is Temperature_data$var. This will the variables saved. E.g. in this TraCE file used in this example temperature is saved as TS for surface temperature.

TS <- ncvar_get(Temperature_data, "TS") # remember to change TS here if you are using a different variable e.g. PRESL / PSL
TS <- TS - 273.15 # here we are transforming the temperature from Kelvin to degreesC

The same thing is true for the time variable, latitude and longitude. Instead of using Temperature_data$var as above, you should use Temperature_data$dim as this will show the dimensions of the data. NOTE: some files will also contain a fourth dimension which is height. In this TraCE temperature file the climate is only simulated at one level (2m above the surface), however other files will contain this. See below an example, however this is not run in this code example.

Time <- ncvar_get(Temperature_data, "time") # units of this are in ka BP (check this by running Temperature_data$dim$time$units)
Time <- Time*-1000 # converting to BP 
lat <- ncvar_get(Temperature_data, "lat")
lon <- ncvar_get(Temperature_data, "lon")
#height <- ncvar_get(Temperature_data, "height") # this is not required for TraCE temperature. 

For this plot we are going to extract the time between 5,000 and 7,000 years before present and create an average. Time is in ka BP and therefore we need to have -5 and -7 to represent the time window.

We are then going to create an average temperature

start_year <- 5000 # this is the year you want it to start (5000 BP)
end_year <- 7000 # end year 
time_indices <- which(Time >= start_year & Time <= end_year)
avg_temperature <- apply(TS[,,time_indices], c(1, 2), mean, na.rm = TRUE)

Now, all of the data required for this plot / extraction is complete. The next step is closing the model file uploaded at the beginning.

NOTE: Because you have closed this file, if you did want to change anything in the script above this (e.g. time period of interest), then you will have to re-run the script!

nc_close(Temperature_data)

Plotting global temperatures

Now we are going to be tailoring the code to plot the specific characteristics we want. In this section of this code we are going to be setting everything up to plot a global view of temperature. We will then go through some different map layouts and changing colours.

Basic global map

This is where we set the dimensions of our map. Here we will be plotting a global map.

We then create a raster called model_raster which includes our dimensions and our model data. Once this is created we then add a coordinate system, here we are using WGS84 and then finally we sometimes need to flip our model data. We do this in the model_raster <- flip(model_raster, direction='y') line. If you use another model you may not have to do this, so just remove if you so wish!

  lon <- seq(-180, 180, length.out = ncol(avg_temperature))
  lat <- seq(-90, 90, length.out = nrow(avg_temperature))
  
  # Create and modify the raster
  model_raster <- raster(t(avg_temperature), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
 
  model_raster <- flip(model_raster, direction='y')  # Flips 
  
  # Convert the raster to points data frame
  raster_points_df <- as.data.frame(rasterToPoints(model_raster))

basic <- ggplot() +
    geom_raster(data = raster_points_df, aes(x = x, y = y, fill = layer)) +
    geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group), fill = NA, color = "black") +
    labs(title = paste("Average Global Temperature from", start_year, "to", end_year, "years BP"),
         x = "Longitude", y = "Latitude", fill = "Temp (°C)") +
    coord_fixed(1.3) +
    theme_minimal()
print(basic)

Changing map colours

The above plot is fine, but the colours do not really resemble temperature. In the plot below we are keeping the same plot, but adding a few lines of code to have a colour gradient.

basic <- ggplot() +
    geom_raster(data = raster_points_df, aes(x = x, y = y, fill = layer)) +
    geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group), fill = NA, color = "black") +
    scale_fill_gradient2(low = "darkblue", mid = "white", high = "red",
                         midpoint = median(raster_points_df$layer, na.rm = TRUE),
                         name = "Temp (°C)", limits = range(raster_points_df$layer, na.rm = TRUE),
                         space = "Lab") +
    labs(title = paste("Average Global Temperature from", start_year, "to", end_year, "years BP"),
         x = "Longitude", y = "Latitude", fill = "Temp (°C)") +
    coord_fixed(1.3) +
    theme_minimal()

print(basic)

Changing map projections

Changing the projection of maps in R can sometimes be a bit challenging and requires a few more lines of code. Here there are two code blocks to prepare the data and world map outline to the correct projection, and then plot

r <- raster(t(avg_temperature), xmn=-180, xmx=180, ymn=-90, ymx=90, crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  r <- flip(r, direction='y')  # This is where we flip vertically again

  robinson_crs <- "+proj=robin +datum=WGS84" # assigning a differnet projection here
  r_robinson <- projectRaster(r, crs=robinson_crs)

  raster_points_df <- as.data.frame(rasterToPoints(r_robinson)) # Converting our data points to a dataframe
  names(raster_points_df) <- c("Longitude", "Latitude", "Temperature") # add column headings

  world_map <- map("world", fill = TRUE, plot = FALSE) # getting world map
  world_sp <- map2SpatialPolygons(world_map, IDs=world_map$names, proj4string=CRS("+proj=longlat +datum=WGS84")) # changing the projections 
  world_sf <- st_as_sf(world_sp)
  robinson_crs <- "+proj=robin +datum=WGS84"
  world_sf <- st_transform(world_sf, crs = robinson_crs)

gg <- ggplot() +
    geom_tile(data = raster_points_df, aes(x = Longitude, y = Latitude, fill = Temperature)) +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red",
                         midpoint = median(raster_points_df$Temperature, na.rm = TRUE), # Adjust midpoint as needed
                         name = "Temperature (°C)") +
    geom_sf(data = world_sf, color = "black", fill = NA) +
    coord_sf() +  
    labs(title = paste("Average Global Temperature from", start_year, "to", end_year, "years BP"),
         x = "Longitude", y = "Latitude", fill = "Temp (°C)") +
    theme_minimal()

print(gg)

Now we have a nice plot of global temperatures. This is not in climate anomalies, but please do change if you so wish, and then in the plotting code above, you may wish to change the area of white to 0 degrees by swapping this line of code: midpoint = median(raster_points_df$Temperature, na.rm = TRUE) to midpoint = 0 .