02: Creating Plots

library(IRdisplay)
display_html('<iframe width="560" height="315" src="https://www.youtube.com/embed/MXr3tp6Q1aA" frameborder="0" allowfullscreen></iframe>')

In this chapter we will look at how to plot data from a NetCDF file. We will look at a 1D data variable as well as some 2D data that we will plot on a map. Let’s load in some libraries that we will use.

if (!requireNamespace("ggplot2", quietly = TRUE)) {
  install.packages("ggplot2")
}
library(ggplot2)
library(RNetCDF)
if (!requireNamespace("maps", quietly = TRUE)) {
  install.packages("maps")
}
library(maps)

1D data

We will start by looking at a depth profile of Chlorophyll A data. If you use these data in a publication, please cite them in your list of references with the following citation:

Anna Vader, Lucie Goraguer, Luke Marsden (2022) Chlorophyll A and phaeopigments Nansen Legacy cruise 2021708 station P4 (NLEG11) 2021-07-18T08:50:42 https://doi.org/10.21335/NMDC-1248407516-P4(NLEG11)

url <- 'https://opendap1.nodc.no/opendap/chemistry/point/cruise/nansen_legacy/2021708/Chlorophyll_A_and_phaeopigments_Nansen_Legacy_cruise_2021708_station_P4_NLEG11_20210718T085042.nc'
data <- open.nc(url)

depth <- var.get.nc(data, 'DEPTH')
chla_total <- var.get.nc(data, 'CHLOROPHYLL_A_TOTAL')

Let’s create a quick line plot.

plot(
  chla_total,
  depth,
  type = "l",
  xlab = "Chlorophyll A",
  ylab = "Depth",
  main = "Line Plot of Depth vs Chla_total",
)

Let’s turn that around so depth is increasing going downwards.

plot(
  chla_total,
  depth,
  type = "l",
  xlab = "Chlorophyll A",
  ylab = "Depth",
  main = "Line Plot of Depth vs Chla_total",
  ylim = c(max(depth), min(depth))
)

We can alternatively make a scatter plot.

plot(
  chla_total,
  depth,
  type = "p",  # Change from "l" to "p" for scatter plot
  xlab = "Chlorophyll A",
  ylab = "Depth",
  main = "Scatter Plot of Depth vs Chla_total",
  ylim = c(max(depth), min(depth))
)

# Enhanced Scatter Plot
plot(
  chla_total,
  depth,
  type = "p",  # Scatter plot
  pch = 19,    # Solid circle points
  col = "blue",# Point color
  cex = 1.5,   # Point size
  xlab = "",
  ylab = "",
  ylim = c(350, 0),
  xlim = c(0, 7),
  las = 1,     # Horizontal y-axis labels
)

# Add a grid for better readability
grid(col = "lightgray", lty = "dotted", lwd = 1)

# Customize axis labels and main title
title(
  main = "Scatter Plot of Depth vs Chlorophyll A",
  col.main = "darkblue",
  font.main = 4,  # Bold and italic main title
  cex.main = 1.2  # Main title size
)
title(
  xlab = "Chlorophyll A (µg/L)",
  ylab = "Depth (m)",
  col.lab = "darkred",
  font.lab = 2,   # Bold axis labels
  cex.lab = 1.2   # Axis labels size
)

# Enhanced Line Plot
plot(
  chla_total,
  depth,
  type = "l",  # Line plot
  lwd = 2,     # Line width
  col = "blue",# Line color
  xlab = "",
  ylab = "",
  ylim = c(350, 0),
  xlim = c(0, 7),
  las = 1     # Horizontal y-axis labels
)

# Add a grid for better readability
grid(col = "lightgray", lty = "dotted", lwd = 1)

# Customize axis labels and main title
title(
  main = "Line Plot of Depth vs Chlorophyll A",
  col.main = "darkblue",
  font.main = 4,  # Bold and italic main title
  cex.main = 1.2  # Main title size
)
title(
  xlab = "Chlorophyll A (µg/L)",
  ylab = "Depth (m)",
  col.lab = "darkred",
  font.lab = 2,   # Bold axis labels
  cex.lab = 1.2   # Axis labels size
)

# Add points on the line plot
points(
  chla_total, 
  depth, 
  pch = 19, 
  col = "red", 
  cex = 1.5
)

Plotting 2D data on a map

We will load in some multidimensional data now. This is a dataset of global surface temperature anomalies from the NOAA Global Surface Temperature Dataset (NOAAGlobalTemp), Version 5.0. We introduced these data in tutorial #01.

https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ncdc:C01585

H.-M. Zhang, B. Huang, J. H. Lawrimore, M. J. Menne, and T. M. Smith (2019): NOAA Global Surface Temperature Dataset (NOAAGlobalTemp), Version 5.0. NOAA National Centers for Environmental Information. doi:10.25921/9qth-2p70 Accessed 2024-04-09.

url <- 'https://www.ncei.noaa.gov/thredds/dodsC/noaa-global-temp-v5/NOAAGlobalTemp_v5.0.0_gridded_s188001_e202212_c20230108T133308.nc'
data <- open.nc(url)

We want to access the anom variable and extract data only for a data of our choice. Let’s repeat the steps we took in tutorial #01 to do this.

latitude <- var.get.nc(data, 'lat')
longitude <- var.get.nc(data, 'lon')

desired_date <- as.Date('2020-01-01')
days_since_1800 <- as.numeric(difftime(desired_date, as.Date('1800-01-01'), units = 'days'))

time <- var.get.nc(data, "time") 
# Finding index of the value
time_index <- which(time == days_since_1800)

anom <- var.get.nc(data, 'anom', start=c(NA, NA, 1, time_index), count=c(NA,NA,1,1))

Now let’s make a very basic plot of the data.

image(longitude, latitude, anom)

That is very basic! Let’s use ggplot.

The longitudes need to be adjusted to between -180° and 180° to align with the coastlines in the map.

Let’s also switch to a white-centred colour palette.

# Convert data to data frame
df <- expand.grid(lon = longitude, lat = latitude)
df$anom <- as.vector(anom)

# Create the plot using ggplot2
ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom)))

Let’s add the coastlines and country borders now.

# Create the plot using ggplot2
ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom))) +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

We need to shift the longitude values in our dataframe to be between -180 and 180 instead of 0 and 360-

# Adjust lon variable
longitude_shifted <- ifelse(longitude > 180, longitude - 360, longitude)

# Convert data to data frame
df <- expand.grid(lon = longitude_shifted, lat = latitude)
df$anom <- as.vector(anom)

# Create the plot using ggplot2
ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom))) +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

We can adjust the aspect ratio

ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom))) +
  coord_fixed(ratio = 1.5) +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

Let’s remove the grey box by using theme_minimal.

ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom))) +
  coord_fixed(ratio = 1.5) +
  theme_minimal() +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

We might also want to change the label on our colour palette to something more understandable. Let’s also move the colour bar to below the axis.

ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom))) +
  coord_fixed(ratio = 1.5) +
  theme_minimal() +
  labs(fill = "Temperature anomaly (°C)") +
  theme(legend.position = "bottom") +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

Finally, let’s add a title.

ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(-10,10)) +
  coord_fixed(ratio = 1.5) +
  theme_minimal() +
  labs(fill = "Temperature anomaly (°C)") +
  theme(legend.position = "bottom") +
  ggtitle(paste("Temperature anomalies ", desired_date)) +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

Maybe we want to zoom in on a certain region.

ggplot(df, aes(x = lon, y = lat, fill = anom)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(min(df$anom), max(df$anom))) +
  coord_fixed(ratio = 1.5) +
  theme_minimal() +
  labs(fill = "Temperature anomaly (°C)") +
  theme(legend.position = "bottom") +
  ggtitle(paste("Temperature anomalies ", desired_date)) +
  xlim(-20, 50) +
  ylim(35, 75) +
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA)

Warning message:
“Removed 2480 rows containing missing values or values outside the scale range (`geom_tile()`).”

Finally, let’s write that as a function so you can easily run it for whichever date you want to. Here’s the whole code including the bits before the function, so you can easily take this with you and run it yourself.

I am going to manually set the limits on the color bar this time so we are using the same range for each plot

library(RNetCDF)
library(ggplot2)

url <- 'https://www.ncei.noaa.gov/thredds/dodsC/noaa-global-temp-v5/NOAAGlobalTemp_v5.0.0_gridded_s188001_e202212_c20230108T133308.nc'
data <- open.nc(url)

plot_temp_anomaly <- function(desired_date, data) {

  latitude <- var.get.nc(data, 'lat')
  longitude <- var.get.nc(data, 'lon')
  longitude_shifted <- ifelse(longitude > 180, longitude - 360, longitude) # Adjust lon variable
  time <- var.get.nc(data, "time")
  
  days_since_1800 <- as.numeric(difftime(desired_date, as.Date('1800-01-01'), units = 'days'))
  
  # Finding index of the value
  time_index <- which(time == days_since_1800)
  
  anom <- var.get.nc(data, 'anom', start=c(NA, NA, 1, time_index), count=c(NA,NA,1,1))
  
  # Convert data to data frame
  df <- expand.grid(lon = longitude_shifted, lat = latitude)
  df$anom <- as.vector(anom)
  
  ggplot(df, aes(x = lon, y = lat, fill = anom)) +
    geom_tile() +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, limits = c(-10,10)) +
    coord_fixed(ratio = 1.5) +
    theme_minimal() +
    labs(fill = "Temperature anomaly (°C)") +
    theme(legend.position = "bottom") +
    ggtitle(paste("Temperature anomalies ", desired_date)) +
    xlim(-180, 180) +
    ylim(-90, 90) +
    geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group),
                 color = "black", fill = NA)
}

# Example usage:
plot_temp_anomaly(as.Date('1960-01-01'), data)
plot_temp_anomaly(as.Date('1990-01-01'), data)
plot_temp_anomaly(as.Date('2020-01-01'), data)

This is only a short demonstration of how easy it can be to plot data out of a CF-NetCDF file. R is great for plotting data and data analysis and there are many good tutorials online, so go and explore!

How to cite this course

If you think this course contributed to the work you are doing, consider citing it in your list of references. Here is a recommended citation:

Marsden, L. (2024, May 31). NetCDF in R - from beginner to pro. Zenodo. https://doi.org/10.5281/zenodo.11400754

And you can navigate to the publication and export the citation in different styles and formats by clicking the icon below.