Climate Visualisations

Followed along with Pat Schloss in his Riffomonas Project.

GISS Surface Temperature Analysis (GISTEMP v4)

The GISS Surface Temperature Analysis ver. 4 (GISTEMP v4) is an estimate of global surface temperature change.

Temperature Index Plot

read_csv("data/GLB.Ts+dSST.csv", skip = 1, na = "***") %>% 
  select(year = Year, t_diff = `J-D`) %>% 
  drop_na() %>% 
  ggplot(aes(x = year, y = t_diff)) +
  geom_line(aes(colour = "1"), size = 0.5, show.legend = FALSE) +
  geom_point(fill = "white", aes(colour = "1"), shape = 21, show.legend = TRUE) +
  geom_smooth(se = FALSE, aes(colour = "2"), size = 0.5, span = 0.2, show.legend = FALSE) + 
  scale_x_continuous(breaks = seq(1880, 2023, 20), expand = c(0,0)) +
  scale_y_continuous(limits = c(-0.5, 1.5), expand = c(0,0)) +
  scale_colour_manual(name = NULL, breaks = c(1,2), values = c("gray", "black"), labels = c("Annual mean", "Loess smoothing"),
                      guide = guide_legend(override.aes = list(shape = 15, size = 5))) +
  labs(x = "", y = "Temperature Anomaly w.r.t 1951-1980 (\u00B0C)",
       title = "Global Annual Mean Surface Air Temperature Change",
       subtitle = "Data source: NASA's Goddard Institute for Space Studies (GISS).\nCredit: NASA/GISS") +
  theme_light() +
  theme(axis.ticks = element_blank(),
        axis.title.y = element_text(size = 9),
        plot.title.position = "plot",
        plot.title = element_text(margin = margin(b = 10), colour = "red", face = "bold"),
        plot.subtitle = element_text(size = 8, margin = margin(b = 10)),
        legend.position = c(0.15, 0.9),
        legend.title = element_text(size = 0),
        legend.key.height = unit(10, "pt"),
        legend.margin = margin(0,0,0,0))

## Rows: 143 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (19): Year, Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

GISTEMP Team, 2022: GISS Surface Temperature Analysis (GISTEMP), version 4. NASA Goddard Institute for Space Studies. Dataset accessed 2022-06-19 at https://data.giss.nasa.gov/gistemp/.

Temperature Bar Plot

t_data <- read_csv("data/GLB.Ts+dSST.csv", skip = 1, na ="***") %>%
  select(year = Year, t_diff = `J-D`) %>% 
  drop_na()

## Rows: 143 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (19): Year, Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

annotation <- t_data %>% 
  arrange(year) %>% 
  slice(1, n()) %>% 
  mutate(t_diff = 0,
         x = year + c(-5,5))

t_data %>% 
  ggplot(aes(x = year, y = t_diff, fill = t_diff)) +
  geom_col(show.legend = FALSE) +
  geom_text(data = annotation, aes(x = x, label = year), colour = "white") +
  geom_text(x = 1880, y = 1, hjust = 0,
            label = glue("Global temperatures have increased by over {max(t_data$t_diff)}\u00B0C since {min(t_data$year)}"),
            colour = "white") +
  scale_fill_stepsn(colours = c("darkblue", "white", "darkred"),
                    values = rescale(c(min(t_data$t_diff), 0, max(t_data$t_diff))),
                    limits = c(min(t_data$t_diff), max(t_data$t_diff)),
                    n.breaks = 9) +
  theme_void() +
  theme(plot.background = element_rect(fill = "black"))

### Temperature Ridgelines Plot

# Download data and unzip
url <- 'https://data.giss.nasa.gov/pub/gistemp/gistemp250_GHCNv4.nc.gz'
download.file(url, destfile = 'gistemp250_GHCNv4.nc.gz')
gunzip('gistemp250_GHCNv4.nc.gz')

nc_data <- nc_open('gistemp250_GHCNv4.nc')

# Save print(nc) dump to a text file
sink('gistemp250_GHCNv4.txt')
print(nc_data)

## File gistemp250_GHCNv4.nc (NC_FORMAT_CLASSIC):
## 
##      2 variables (excluding dimension variables):
##         int time_bnds[nv,time]   
##         short tempanomaly[lon,lat,time]   
##             long_name: Surface temperature anomaly
##             units: K
##             scale_factor: 0.00999999977648258
##             cell_methods: time: mean
##             _FillValue: 32767
## 
##      4 dimensions:
##         lat  Size:90 
##             standard_name: latitude
##             long_name: Latitude
##             units: degrees_north
##         lon  Size:180 
##             standard_name: longitude
##             long_name: Longitude
##             units: degrees_east
##         time  Size:1715 
##             long_name: time
##             units: days since 1800-01-01 00:00:00
##             bounds: time_bnds
##         nv  Size:2 (no dimvar)
## 
##     5 global attributes:
##         title: GISTEMP Surface Temperature Analysis
##         institution: NASA Goddard Institute for Space Studies
##         source: http://data.giss.nasa.gov/gistemp/
##         Conventions: CF-1.6
##         history: Created 2022-12-12 11:34:41 by SBBX_to_nc 2.0 - ILAND=250,  IOCEAN=none,     Base: 1951-1980

sink()

lon <- ncvar_get(nc_data, 'lon')
lat <- ncvar_get(nc_data,'lat', verbose = FALSE)
t <- ncvar_get(nc_data, 'time')

# Store data in 3D array
t_anomaly.array <- ncvar_get(nc_data, 'tempanomaly')
dim(t_anomaly.array)

## [1]  180   90 1715

fillvalue <- ncatt_get(nc_data, 'tempanomaly', '_FillValue')
t_anomaly.array[t_anomaly.array == fillvalue$value] <- NA

# Wrangle data
# Convert 3D array into tibble
t_data <- as.data.table(t_anomaly.array) %>%
  as_tibble() %>% 
  select(longitude = V1, latitude = V2, time = V3, t_diff = value) %>% 
  mutate(longitude = lon[longitude],
         latitude = lat[latitude],
         time = t[time] + as.Date('1800-01-01'),
         year = year(time)) %>% 
  group_by(year, longitude, latitude) %>% 
  summarise(t_diff = mean(t_diff), .groups = 'drop') %>% 
  filter(year >= 1950 & year < 2022) %>% 
  group_by(year) %>% 
  mutate(t_ave = mean(t_diff))

# Plot
t_data %>% 
  ggplot(aes(x = t_diff, y = factor(year, levels = seq(2021, 1950, -1)), 
             fill = t_ave)) +
  geom_density_ridges(bandwidth = 0.2, scale = 3.5, size = 0.2) +
  scale_fill_gradient2(low = 'darkblue', mid = 'white', high = 'darkred',
                       midpoint = 0) +
  coord_cartesian(xlim = c(-4.5, 4.5)) +
  scale_x_continuous(breaks = seq(-4, 4, 2)) +
  scale_y_discrete(breaks = seq(1950, 2021, 10)) +
  labs(x = 'Temperature anomaly (\u00B0 C)',
       y = '',
       title = 'Land Temperature Anomalies') +
  theme(panel.background = element_rect(fill = 'black'),
        plot.background = element_rect(fill = 'black'),
        panel.grid = element_blank(),
        text = element_text(colour = 'white'),
        axis.text = element_text(colour = 'white'),
        axis.ticks = element_line(colour = 'white'),
        axis.ticks.y = element_blank(),
        axis.line.x = element_line(colour = 'white'),
        axis.line.y = element_blank(),
        legend.position = 'none')

# Unlink ncdata
nc_close(nc_data)
unlink('gistemp250_GHCNv4.nc')
unlink('gistemp250_GHCNv4.txt')

Temperature Anomalies Map

# Download data and unzip
url <- 'https://data.giss.nasa.gov/pub/gistemp/gistemp250_GHCNv4.nc.gz'
download.file(url, destfile = 'gistemp250_GHCNv4.nc.gz')
gunzip('gistemp250_GHCNv4.nc.gz')

nc_data <- nc_open('gistemp250_GHCNv4.nc')

# Save print(nc) dump to a text file
sink('gistemp250_GHCNv4.txt')
print(nc_data)

## File gistemp250_GHCNv4.nc (NC_FORMAT_CLASSIC):
## 
##      2 variables (excluding dimension variables):
##         int time_bnds[nv,time]   
##         short tempanomaly[lon,lat,time]   
##             long_name: Surface temperature anomaly
##             units: K
##             scale_factor: 0.00999999977648258
##             cell_methods: time: mean
##             _FillValue: 32767
## 
##      4 dimensions:
##         lat  Size:90 
##             standard_name: latitude
##             long_name: Latitude
##             units: degrees_north
##         lon  Size:180 
##             standard_name: longitude
##             long_name: Longitude
##             units: degrees_east
##         time  Size:1715 
##             long_name: time
##             units: days since 1800-01-01 00:00:00
##             bounds: time_bnds
##         nv  Size:2 (no dimvar)
## 
##     5 global attributes:
##         title: GISTEMP Surface Temperature Analysis
##         institution: NASA Goddard Institute for Space Studies
##         source: http://data.giss.nasa.gov/gistemp/
##         Conventions: CF-1.6
##         history: Created 2022-12-12 11:34:41 by SBBX_to_nc 2.0 - ILAND=250,  IOCEAN=none,     Base: 1951-1980

sink()

lon <- ncvar_get(nc_data, 'lon')
lat <- ncvar_get(nc_data,'lat', verbose = FALSE)
t <- ncvar_get(nc_data, 'time')

# Store data in 3D array
t_anomaly.array <- ncvar_get(nc_data, 'tempanomaly')
dim(t_anomaly.array)

## [1]  180   90 1715

fillvalue <- ncatt_get(nc_data, 'tempanomaly', '_FillValue')
t_anomaly.array[t_anomaly.array == fillvalue$value] <- NA

# Wrangle data
# Convert 3D array into tibble
t_data <- as.data.table(t_anomaly.array) %>%
  as_tibble() %>% 
  select(longitude = V1, latitude = V2, time = V3, t_diff = value) %>% 
  mutate(longitude = lon[longitude],
         latitude = lat[latitude],
         time = t[time] + as.Date('1800-01-01'),
         year = year(time)) %>% 
  group_by(year, longitude, latitude) %>% 
  summarise(t_diff = mean(t_diff), .groups = 'drop') %>% 
  filter(year >= 1950 & year < 2022) %>% 
  mutate(decade = 10 * floor(year / 10),
         single = year %% 10)
  
# Plot temperature raster maps
t_data %>%
  mutate(t_diff = case_when(t_diff < -4 ~ -4,
                            t_diff > 4 ~ 4,
                            TRUE ~ t_diff)) %>%
  ggplot(aes(x = longitude, y = latitude, fill = t_diff)) +
  geom_raster() +
  facet_grid(decade~single, switch = 'y') +
  scale_fill_gradient2(name = 'Temperature anomaly (\u00B0 C)',
                       low = 'darkblue', mid = 'white', high = 'darkred',
                       midpoint = 0,
                       limits = c(-5, 5),
                       breaks = c(-4, -2, 0, 2, 4)) +
  coord_fixed(expand = FALSE) +
  labs(x = '', y = '',
       title = 'Global Temperature Anomalies (1950 - 2021)') +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_rect(fill = 'black'),
        plot.background = element_rect(fill = 'black'),
        panel.grid = element_blank(),
        strip.text.x = element_blank(),
        strip.text.y.left = element_text(angle = 0, colour = 'white'),
        strip.background = element_blank(),
        title = element_text(colour = 'white', face = 'bold'),
        legend.position = c(0.75, 0.04),
        legend.direction = 'horizontal',
        legend.title = element_text(colour = 'white', size = 6),
        legend.text = element_text(colour = 'white', size = 5),
        legend.background = element_rect(fill = 'black')) +
  guides(fill = guide_colorbar(title.position = 'top',
                               title.hjust = 0.5))

# Unlink ncdata
nc_close(nc_data)
unlink('gistemp250_GHCNv4.nc')
unlink('gistemp250_GHCNv4.txt')

Temperature Rug Plot

# Load data
zone_data <- read_csv("https://data.giss.nasa.gov/gistemp/tabledata_v4/ZonAnn.Ts+dSST.csv")

## Rows: 142 Columns: 15
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (15): Year, Glob, NHem, SHem, 24N-90N, 24S-24N, 90S-24S, 64N-90N, 44N-64...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Wrangle data
bands <- rev(c("64N-90N", "44N-64N", "24N-44N", "EQU-24N", "24S-EQU", "44S-24S", "64S-44S", "90S-64S"))

zone_data <- read_csv("https://data.giss.nasa.gov/gistemp/tabledata_v4/ZonAnn.Ts+dSST.csv") %>%
  select(year = Year, all_of(bands)) %>%
  pivot_longer(-year, names_to = "zone", values_to = "t_diff") %>%
  mutate(zone = factor(zone, levels = bands),
         zone_position = as.numeric(zone))

## Rows: 142 Columns: 15
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (15): Year, Glob, NHem, SHem, 24N-90N, 24S-24N, 90S-24S, 64N-90N, 44N-64...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

current_year <- zone_data %>% filter(year == 2021)

# Plot
zone_data %>%
  ggplot(aes(x = t_diff, xend = t_diff,
             y = zone_position - 0.25, yend = zone_position + 0.25)) +
  geom_segment(colour = "white", alpha = 0.25) +
  geom_segment(data = current_year,
               aes(colour = t_diff), size = 2, lineend = "round") +
  scale_y_continuous(breaks = 1:8,
                     labels = bands) +
  scale_x_continuous(breaks = seq(-3, 4, 1),
                     labels = seq(-3, 4, 1),
                     limits = c(-3, 4)) +
  scale_colour_gradient2(low = "darkblue", mid = "white", high = "darkred",
                         midpoint = 0, guide = "none") +
  labs(x = "Temperature anomaly (\u00B0C)", y = "",
       title = "Temperature Anomaly by Latitude (1880 - 2021)",
       subtitle = "Bars for 2021 are coloured by the size of the anomaly") +
  theme(plot.background = element_rect(fill = "black", colour = "black"),
        panel.background = element_rect(fill = "black"),
        plot.title = element_text(colour = "white", face = "bold"),
        plot.subtitle = element_text(colour = "gray", size = 9),
        plot.title.position = "plot",
        axis.text = element_text(colour = "white"),
        axis.title = element_text(colour = "white"),
        panel.grid.major.x = element_line(colour = "gray", size = 0.25),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_line(colour = "white"),
        axis.ticks = element_blank())