Make money from data - CYBAEA.NET - Area Plots with Intensity Coloring

I am not sure apeescape’s ggplot2 area plot with intensity colouring (dead link) is really the best way of presenting the information, but it had me intrigued enough to replicate it using base R graphics.

The key technique is to draw a gradient line which R does not support natively so we have to roll our own code for that. Unfortunately, lines(..., type="l") does not recycle the colour col= argument, so we end up with rather more loops than I thought would be necessary.

(The answer is not to use lines(..., type="h") which, confusingly, does recycle the colour col= argument. This one had me for a while, but the type=h lines always start from zero so you do not get the gradient feature.)

We also get a nice opportunity to use the under-appreciated read.fwf function.

Code

##!/usr/bin/Rscript
## nino.R - another version of http://bit.ly/9P9Gh1
## Copyright © 2010 Allan Engelhardt (http://www.cybaea.net/)

## Get the data from the NOAA server
nino <- read.fwf(
  "https://www.cpc.ncep.noaa.gov/data/indices/wksst.for",
  widths = c(-1, 9, rep(c(-5, 4, 4), 4)),
  skip = 4,
  col.names = c(
    "Week",
    paste(
      rep(
        c("Nino12", "Nino3", "Nino34", "Nino4"), rep(2, 4)
      ),
      c("SST", "SSTA"),
      sep = "."
    )
  )
)

## Make the date column something useful
nino$Week <- as.Date(nino$Week, format = "%d%b%Y")

## Make colour gradients
ncol <- 50
grad.neg <- hsv(4 / 6, seq(0, 1, length.out = ncol), 1) # Blue gradient
grad.pos <- hsv(0, seq(0, 1, length.out = ncol), 1) # Red gradient

## Make plot
{ ## <- This block to work in R Markdown / Quarto
  plot(
    Nino34.SSTA ~ Week,
    data = nino,
    type = "n",
    main = "Nino34",
    xlab = "Date",
    ylab = "SSTA",
    axes = FALSE
  )
  do.call(
    function(...) {
      rect(..., col = "gray85", border = NA)
    },
    as.list(par("usr")[c(1, 3, 2, 4)])
  )

  y <- nino$Nino34.SSTA # The values we will plot
  x <- nino$Week

  axis.Date(1,
    x = x,
    tck = 1,
    col = "white"
  )
  axis(2, tck = 1, col = "white")
  box()

  idx <- integer(NROW(nino))
  idx[y >= 0] <-
    1 + round(y[y >= 0] * (ncol - 1) / max(y[y >= 0]), 0)
  idx[y < 0] <-
    1 + round(-y[y < 0] * (ncol - 1) / max(-y[y < 0]), 0)

  draw.gradient <- function(x, ys, cols) {
    xs <- rep(x, 2)
    for (i in seq(1, length(ys) - 1)) {
      plot.xy(list(x = xs, y = c(ys[i], ys[i + 1])),
        type = "l", col =
          cols[i]
      )
    }
  }

  for (i in 1:length(x)) {
    ys <- seq(0, y[i], length.out = idx[i] + 1)
    cols <- (if (y[i] >= 0) {
      grad.pos
    } else {
      grad.neg
    })
    draw.gradient(x[i], ys, cols)
  }
}

The result is a decent gradient:

I deliberately omitted the scale legend on the right hand side following Allan’s First Law of Happy Graphics: Thou shall not present the same information twice.

For less dense information, you should increase the line width. That is left to the reader. (Hint: it is hard to get just right in base graphics, but lwd <- ceiling(par("pin")[1] / dev.size("in")[1] * dev.size("px")[1] / length(x)) could be a starting point for an approximation. We really need gradient-filled polygons in base R.)