Add color gradient to ridgelines according to height

I want to color my ridgeline plot with gradient fill colors depending on the height of the area instead of depending on the X axis. It would be somthing like this:

library(ggplot2)
library(ggridges)   
ggplot(lincoln_weather, aes(x = `Mean Temperature [F]`, y = `Month`, fill = stat(x))) +
      geom_density_ridges_gradient(scale = 3, size = 0.3, rel_min_height = 0.01) +
      scale_fill_viridis_c(name = "Temp. [F]", option = "C") +
      labs(title = 'Temperatures in Lincoln NE')

But here instead of fill the gradient according to the x axis I'm looking to use the same vertical gradient for all the curves so the highest the curve, the darkest the color. So the peaks would have a darker color while near the baseline would be lighter.

PD:my plot I use geom_ridgeline() instead of geom_density_ridges_gradient() but I think this example is more clear to illustrate the problem. And also have negative values.

PD2:I know there is a similar question "solved" here. But is obsolete, since if you check the github from the suggested package there is an issue stating that the function needed doesn't work because is based on another deprecated function.

An option often overlooked by people is that you can pretty much draw anything in ggplot2 as long as you can express it in polygons. The downside is that it is a lot more work.

The following approach mimics my answer given here but adapted a little bit to work more consistently for multiple densities. It abandons the {ggridges} approach completely.

The function below can be used to slice up an arbitrary polygon along the y-position.

library(ggplot2)
library(ggridges)
library(polyclip)

fade_polygon <- function(x, y, yseq = seq(min(y), max(y), length.out = 100)) {
  poly <- data.frame(x = x, y = y)
  
  # Create bounding-box edges
  xlim <- range(poly$x) + c(-1, 1)
  
  # Pair y-edges
  grad <- cbind(head(yseq, -1), tail(yseq, -1))
  # Add vertical ID
  grad <- cbind(grad, seq_len(nrow(grad)))
  
  # Slice up the polygon
  grad <- apply(grad, 1, function(range) {
    # Create bounding box
    bbox <- data.frame(x = c(xlim, rev(xlim)),
                       y = c(range[1], range[1:2], range[2]))
    
    # Do actual slicing
    slice <- polyclip::polyclip(poly, bbox)
    
    # Format as data.frame
    for (i in seq_along(slice)) {
      slice[[i]] <- data.frame(
        x = slice[[i]]$x,
        y = slice[[i]]$y,
        value = range[3],
        id = c(1, rep(0, length(slice[[i]]$x) - 1))
      )
    }
    slice <- do.call(rbind, slice)
  })
  # Combine slices
  grad <- do.call(rbind, grad)
  # Create IDs
  grad$id <- cumsum(grad$id)
  return(grad)
}

Next, we need to calculate densities manually for every month and apply the function above for each one of those densities.

# Split by month and calculate densities
densities <- split(lincoln_weather, lincoln_weather$Month)
densities <- lapply(densities, function(df) {
  dens <- density(df$`Mean Temperature [F]`)
  data.frame(x = dens$x, y = dens$y)
})

# Extract x/y positions
x <- lapply(densities, `[[`, "x")
y <- lapply(densities, `[[`, "y")

# Make sequence to max density
ymax <- max(unlist(y))
yseq <- seq(0, ymax, length.out = 100) # 100 can be any large enough number

# Apply function to all densities
polygons <- mapply(fade_polygon, x = x, y = y, yseq = list(yseq),
                   SIMPLIFY = FALSE)

Next, we need to add the information about Months back into the data.

# Count number of observations in each of the polygons
rows <- vapply(polygons, nrow, integer(1))
# Combine all of the polygons
polygons <- do.call(rbind, polygons)
# Assign month information
polygons$month_id  <- rep(seq_along(rows), rows)

Lastly we plot these polygons with vanilla ggplot2. The (y / ymax) * scale does a similar scaling to what ggridges does and adding the month_id offsets each month from oneanother.

scale <- 3
ggplot(polygons, aes(x, (y / ymax) * scale + month_id, 
                     fill = value, group = interaction(month_id, id))) +
  geom_polygon(aes(colour = after_scale(fill)), size = 0.3) +
  scale_y_continuous(
    name = "Month",
    breaks = seq_along(rows),
    labels = names(rows)
  ) +
  scale_fill_viridis_c()

^{Created on 2021-09-12 by the reprex package (v2.0.1)}