Land Temperature Anomaly Distribution

Introduction

Using NASA data, we create an animated plot to show the distribution of the “Land Temperature Anomaly” since the 1960s. The original version is here and in this Medium post. I simply got the original code and translated it into R. The graphics are not as smooth as in the original Python version, but it’s a starting point and you can keep working on it.

Libraries

In this section with install/call all the necessary libraries to process the data

# Load necessary libraries
library(ncdf4)
library(httr)
library(R.utils)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(magick)

The echo: false option disables the printing of code (only output is displayed).

Loading the data

# --- Data Loading ---
url <- "https://data.giss.nasa.gov/pub/gistemp/gistemp250_GHCNv4.nc.gz"
temp_gz <- "gistemp.nc.gz"
temp_nc <- "gistemp.nc"

# Download and decompress
GET(url, write_disk(temp_gz, overwrite = TRUE))

Response [https://data.giss.nasa.gov/pub/gistemp/gistemp250_GHCNv4.nc.gz]
  Date: 2026-02-23 18:15
  Status: 200
  Content-Type: application/x-gzip
  Size: 11.4 MB
<ON DISK>  gistemp.nc.gz

gunzip(temp_gz, destname = temp_nc, remove = TRUE, overwrite = TRUE)

# Open NetCDF
nc_data <- nc_open(temp_nc)
v <- ncvar_get(nc_data, "tempanomaly")
lats <- ncvar_get(nc_data, "lat")
times <- ncvar_get(nc_data, "time") # Usually months since 1800-01-01

Distribution Calculation

# --- Distribution Calculation ---
tOffset <- 82 
timeStep <- 1 
timeBins <- 1 
plotRange <- c(-5, 5)
xx <- seq(plotRange[1], plotRange[2], length.out = 2048)

# Latitude weights (Area weighting)
lat_weights_vec <- cos(lats * pi / 180)

Function to created the base PDF files

makePDFs <- function() {
  nTimeSteps <- floor((dim(v)[3] - 12 * tOffset) / (12 * timeStep))
  pdf_list <- list()
  dates <- c()
  
  for (i in 0:(nTimeSteps - 1)) {
    idxMin <- 12 * (i * timeStep + tOffset) + 1
    idxMax <- 12 * (i * timeStep + timeBins + tOffset)
    
    # Extract slice and weights
    slice <- v[, , idxMin:idxMax]
    # Replicate weights to match slice dimensions
    w_slice <- array(lat_weights_vec, dim = dim(slice)[1:2])
    w_slice <- array(w_slice, dim = dim(slice))
    
    # Flatten and remove NAs (equivalent to .compressed() in Python)
    valid_indices <- !is.na(slice)
    data_points <- slice[valid_indices]
    weights_points <- w_slice[valid_indices]
    
    # Calculate weighted KDE
    # R's density() doesn't support weights as easily as scipy, 
    # so we use the 'weights' argument in density()
    d <- density(data_points, weights = weights_points/sum(weights_points), 
                 from = plotRange[1], to = plotRange[2], n = 2048)
    
    pdf_list[[i + 1]] <- d$y
    dates[i + 1] <- 1880 + (idxMin - 1) / 12
  }
  
  return(list(pdfs = pdf_list, dates = dates))
}



results <- makePDFs()
PDFs <- results$pdfs
startDates <- results$dates

Set colors for animated chart

# --- Setup Colormap ---
# Note: In R, we usually define the gradient manually or use a palette
# Here we approximate the GISS colormap (Blue -> White -> Red)

#giss_colors <- colorRampPalette(c("#0000FF", "#FFFFFF", "#FF0000"))(254)

giss_colors <- colorRampPalette(c("dodgerblue4", "steelblue1", "#FFFFFF", "orange", "indianred"))(512)

Plotting Function

# --- Plotting Function ---
GISSfig <- function(plotTime, startDates, PDFs, metric = TRUE) {
  # Interpolation logic
  nPlot <- floor((plotTime - startDates[1]) / timeStep) + 1
  interpFrac <- (plotTime - startDates[1]) %% timeStep
  
  # Linear interpolation between years
  current_pdf <- PDFs[[nPlot]] * (1 - interpFrac) + PDFs[[nPlot + 1]] * interpFrac
  
  df <- data.frame(x = xx, y = current_pdf)
  
  # Conversion for Fahrenheit
  if (!metric) {
    df$x <- df$x * (9/5)
    # Note: KDE density values must be scaled if the axis is transformed 
    # to maintain area=1, but for visualization we often keep the height.
  }
  
  p <- ggplot(df, aes(x = x, y = y)) +
    # Fill with gradient
    geom_area(aes(fill = x)) +
    scale_fill_gradientn(colors = giss_colors, limits = if(metric) c(-2, 2) else c(-3.6, 3.6), oob = squish) +
    # Main line
    geom_line(color = "white", linewidth = 1) +
    # Styling
    theme_dark() +
    theme(
      panel.background = element_rect(fill = "black"),
      plot.background = element_rect(fill = "black"),
      panel.grid.major = element_line(color = "gray30"),
      panel.grid.minor = element_blank(),
      axis.text.y = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(size = 22),
      axis.title.x = element_text(size = 22),
      plot.title = element_text(size = 30, face = "bold"),
      text = element_text(color = "white")
    ) +
    labs(
      title = "Land Temperature Anomaly Distribution",
      x = if(metric) "Temperature Anomaly (°C)" else "Temperature Anomaly (°F)"
    ) +
    annotate("text", x = max(df$x)*0.8, y = 0.5, label = as.character(floor(plotTime)), 
             color = "white", size = 16) +
    coord_cartesian(xlim = if(metric) c(-5, 5) else c(-9, 9), ylim = c(0, 0.6))
  
  return(p)
}

Animation Loop

# --- Animation Loop (Example for first 20 frames) --- Now is 120
# In R, for full animation, you'd typically use the 'gganimate' package 
# or loop and save with 'ggsave'.
dir.create("frames_R", showWarnings = FALSE)

figTimes <- seq(1962, 2022.99, length.out = 120) # Reduced for example
for (i in seq_along(figTimes)) {
  p <- GISSfig(figTimes[i], startDates, PDFs, metric = FALSE)
  ggsave(sprintf("frames_R/GISS_%03d.png", i), plot = p, width = 16, height = 9, dpi = 120)
}

Create gif File

## Animation

# List files in the directory, ensure they are sorted correctly
# Use full.names = TRUE to get the complete file paths
img_files <- list.files(path = "frames_R", pattern = "*.png", full.names = TRUE)

# Read all images into a single magick image object
img_list <- lapply(img_files, image_read)
img_joined <- image_join(img_list)


# Animate at 10 frames per second (fps)
img_animated <- image_animate(img_joined, fps = 20)

# Save to disk as a GIF
image_write(image = img_animated, path = "frames_R/my_animation.gif")

Display

knitr::include_graphics("frames_R/my_animation.gif")