# Area Plots with Intensity Coloring

R
Published

13 July 2010

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

## Get the data from the NOAA server
"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")

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) {
} else {
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") / dev.size("in") * dev.size("px") / length(x))` could be a starting point for an approximation. We really need gradient-filled polygons in base R.)