ggplot code for drawing a spectrum with frequency as x-axis, visible bands coloured

rainbow_spectrum.R

library(readxl)
library(dplyr)
library(ggplot2)
library(ragg)
## https://www.sthda.com/english/wiki/ggplot2-themes-and-background-colors-the-3-elements#google_vignette
theme_set(theme_classic(base_size = 16) +
            theme(panel.background = element_rect(fill = "lightgray")))


library(scales) # trans_new() is in the scales library
transform_reciprocal <-
  function() trans_new("reciprocal", function(x) 1/x, function(x) 1/x,
                       domain = c(0, Inf))

## https://www.r-bloggers.com/2012/08/custom-axis-transformations-in-ggplot2/
## https://academo.org/demos/wavelength-to-colour-relationship/
## https://www.johndcook.com/wavelength_to_RGB.html

## FIXME: this could be somewhat vectorized (roughly the same operations
##  for every bin/wavelength range)
nmToRGB <- function(wavelength, Gamma = 0.8, IntensityMax = 255) {
  res <- if(380 <= wavelength && (wavelength<440)) {
    c(red = -(wavelength - 440) / (440 - 380), 0, 1)
  } else if ((wavelength >= 440) && (wavelength<490)){
    c(red = 0, green = (wavelength - 440) / (490 - 440), 
      blue = 1.0)
  } else if((wavelength >= 490) && (wavelength<510)){
    c(red = 0.0,
      green = 1.0,
      blue = -(wavelength - 510) / (510 - 490))
  } else if((wavelength >= 510) && (wavelength<580)){
    c(red = (wavelength - 510) / (580 - 510),
      green = 1.0,
      blue = 0.0)
  } else if((wavelength >= 580) && (wavelength<645)){
    c(red = 1.0, 
      green = -(wavelength - 645) / (645 - 580),
      blue = 0.0)
  } else if((wavelength >= 645) && (wavelength<781)){
    c(red = 1.0,
      green = 0.0,
      blue = 0.0)
  } else {
    c(red = 0.0, green = 0.0, blue = 0.0)
  }
  ## Let the intensity fall off near the vision limits
  if((wavelength >= 380) && (wavelength<420)){
    factor = 0.3 + 0.7*(wavelength - 380) / (420 - 380)
  }else if((wavelength >= 420) && (wavelength<701)){
    factor = 1.0
  } else if((wavelength >= 701) && (wavelength<781)){
    factor = 0.3 + 0.7*(780 - wavelength) / (780 - 700);
  }else{
    factor = 0.0;
  }
  res <- ifelse(res == 0, 0,
                round(IntensityMax * (res * factor)^Gamma)
                )
  res
}

nmToHex <- function(x) {
  vapply(x,
         \(x) do.call(rgb, c(as.list(nmToRGB(x)), maxColorValue = 255)),
         FUN.VALUE = character(1))
}

nmToHex(550:560)

## Download the solar spectrum data
solar_data <- "solar.xls"
if (!file.exists(solar_data)) {
  url <- "https://www.pveducation.org/sites/default/files/PVCDROM/Appendices/AM0AM1_5.xls"
  download.file(url, dest = solar_data)
}

## process
dd <- read_excel(solar_data, col_names = FALSE, skip  = 2, sheet = "Spectra") |>
  setNames(c("Wavelength1",
             paste0("AM1.5_", c("extraterr", "global_tilt", "direct_circ")),
             "blank",
             "Wavelength2", "AM0")) |>
  select(wavelen = Wavelength2, pwr = AM0) |>
  filter(!is.na(wavelen),
                wavelen < 4000) |>
  mutate(colour = nmToHex(wavelen),
                mid_lwr = c(119, (wavelen[-1] + wavelen[-length(wavelen)]) / 2),
         mid_upr = lead(mid_lwr),
         next_wavelen = lead(wavelen))


gg0 <- ggplot(dd, aes(y = pwr)) +
  scale_fill_identity() +
  scale_y_continuous(expand = c(0,0))

## plot by wavelength
gg_WL <- gg0 +
  scale_x_continuous(limits=c(0, 4000), expand = c(0,0)) +
  geom_rect(aes(xmin = mid_lwr, xmax = mid_upr,
                ymin = 0, ymax = pwr, fill = colour))

print(gg_WL)

## plot by frequency
gg_freq <- gg0 +
  ## scale_x_continuous(limits=c(0, 4000), expand = c(0,0)) +
  geom_rect(aes(xmin = 1/mid_lwr, xmax = 1/mid_upr,
                ymin = 0, ymax = pwr, fill = colour)) +
    scale_x_continuous(limits=c(0, 0.005), expand = c(0,0))

print(gg_freq)


## get rid of moire pattern
dd_trunc <- dd |> filter(colour != "#000000")

## filling in all bars gives moire pattern in black region
## using geom_ribbon gives some slight black edges (?because of linear interpolation in geom_ribbon?
##   figure out how to use geom_trapezoid/geom_poly for linear interpolation with bars?)

gg_freq2 <- gg0 +
  geom_ribbon(fill = "black", aes(x = wavelen, ymin=0, ymax = pwr)) +
  scale_x_continuous(expand = c(0,0),
                     transform = "reciprocal",
                     limits =c(4000, 200),
                     breaks = c(100, seq(300, 700, by = 100), 1000)) +
  geom_rect(data = dd_trunc,
            aes(xmin = mid_lwr, xmax = mid_upr,
                ymin = 0, ymax = pwr, fill = colour)) +
  labs(x = "Frequency\n(scale in wavelength, nm)",
       y = expression(atop("Power spectral density",W %.% m^{-2} %.% "nm"^{-1})))

## gg_WL + coord_trans(x = "reciprocal")
## gg_WL + scale_x_continuous(trans = "reciprocal")

ggsave(plot = gg_freq2, "rainbow_frequency.pdf", width = 8, height = 6)
ggsave(plot = gg_freq2, "rainbow_frequency.png")
ggsave(plot = gg_freq2, "rainbow_frequency.png", device = ragg::agg_png, units = "in", width = 8, height = 6)

## how do we get rid of moiré/other effects?

## ideally would like to plot these as a trapezoid: will be annoying to do the necessary data
## manipulation (each row will have to get expanded into a set of 4 rows ...)

## https://ggplot2.tidyverse.org/articles/extending-ggplot2.html?q=new%20geom#creating-a-new-geom

## expand/reshape data into four lines per lines: vertices of polygon?