ggplot_bivariate_diy.R

library(tidyverse)
library(prismatic)

d <- expand.grid(x = 1:25, y = 1:25) 

bi_palette <- function(x, y, palette = viridis::viridis(5), dir_x = 1, dir_y = 1, y_fun = prismatic::clr_desaturate, limits_x = NULL, limits_y = NULL){
  stopifnot(length(x) == length(y))
  if(is.null(limits_x)){ limits_x <- range(x * dir_x, na.rm = TRUE) }
  if(is.null(limits_y)){ limits_y <- range(y * dir_y, na.rm = TRUE) }
  clr_1 <- scales::colour_ramp(palette)(scales::rescale(x * dir_x, from = limits_x))
  map2_chr(.x = clr_1, .y = scales::rescale(y * dir_y, from = limits_y), .f = y_fun)
}

bi_palette_2 <- function(x, y, palette1 = RColorBrewer::brewer.pal(5,"Blues"), palette2 = RColorBrewer::brewer.pal(5, "Oranges"), dir_x = 1, dir_y = 1, limits_x = NULL, limits_y = NULL){
  stopifnot(length(x) == length(y))
  if(is.null(limits_x)){ limits_x <- range(x * dir_x, na.rm = TRUE) }
  if(is.null(limits_y)){ limits_y <- range(y * dir_y, na.rm = TRUE) }
  clr_1 <- scales::colour_ramp(palette1)(scales::rescale(x * dir_x, from = limits_x))
  clr_2 <- scales::colour_ramp(palette2)(scales::rescale(y * dir_y, from = limits_y))
  map2_chr(.x = clr_1, .y = clr_2, .f = prismatic::clr_mix)
}

bi_palette(x = 1:20, y = 26:7) %>%
  color()

ggplot(d, aes(x, y)) +
  geom_tile(aes(fill1 = x,
                fill2 = y,
                fill = after_stat(bi_palette_2(fill1, fill2,
                                               palette1 = RColorBrewer::brewer.pal(5, "Reds"),
                                               palette2 = RColorBrewer::brewer.pal(5,"PuOr")
                )
                ),
                color = after_scale(clr_darken(fill))),
            width = .93, height = .93, size = .5) +
  scale_fill_identity()+
  coord_equal() +
  theme_void() +
  theme(legend.position = "none")


clr <- rcartocolor::carto_pal(7, "Geyser")
expand.grid(x = 1:15, y = 1:15) %>% 
  ggplot(aes(x, y)) +
  geom_tile(aes(fill1 = x,
                fill2 = y,
                fill = after_stat(bi_palette(fill1, fill2,palette = clr, y_fun = clr_desaturate)),
                color = after_scale(clr_darken(fill))),
            width = .93, height = .93, size = .5) +
  scale_fill_identity()+
  coord_equal() +
  theme_void() +
  theme(legend.position = "none")

legend_bivariate <- function(x = c(1, 15),y = c(1, 15), x_title = "x", y_title = "y",dir_x = 1,dir_y = 1, pal = clr, fun, n_bins = 15){
  data <- expand.grid(x_d = 1:n_bins,
                      y_d = 1:n_bins)
  
  x_r <- range(x, na.rm = TRUE)
  y_r <- range(y, na.rm = TRUE)
  x_lab <- x_r %>% scales::breaks_pretty()(3)
  y_lab <- y_r %>% scales::breaks_pretty()(3)
  x_br <- scales::rescale(x_lab, from = x_r, to = c(1,n_bins))
  y_br <- scales::rescale(y_lab, from = y_r, to = c(1,n_bins))
  x_f <- dplyr::between(x_br, left = 1, right = n_bins)
  y_f <- dplyr::between(y_br, left = 1, right = n_bins)

    data %>% 
    ggplot(aes(x = x_d, y = y_d)) +
    geom_tile(aes(fill1 = x_d,
                  fill2 = y_d,
                  fill = after_stat(bi_palette(fill1, fill2,
                                               palette = pal,
                                               y_fun = clr_desaturate,dir_x = dir_x,dir_y = dir_y)))) +
    scale_fill_identity() +
    scale_x_continuous(name = x_title, breaks = x_br, labels = x_lab) +
    scale_y_continuous(name = y_title, breaks = y_br, labels = y_lab) +
    coord_equal(expand = 0) +
    # theme_void() +
    theme_minimal() +
    theme(axis.title = element_text(margin = margin(5,5,5,5)),
          axis.line = element_line(),
          axis.ticks = element_line(),
          panel.grid = element_blank(),
          legend.position = "none")
}
# use case
library(sf)
library(patchwork)
us <- read_sf("~/work/geo_store/USA/usa_states_albers_revised.gpkg")


ggplot() +
  geom_sf(data = us ,
          aes(fill1 = median_income,
              fill2 = log10(popdens),
              fill = after_stat(bi_palette(fill1, fill2,dir_y = -1, palette = clr)),
              color = after_scale(clr_darken(fill)))) +
  scale_fill_identity() +
  theme_minimal() +
  coord_sf(expand = 0,
           crs = 5070) +
  legend_bivariate(x = us$median_income,
                   y = log10(us$popdens),
                   dir_y = -1,
                   pal = clr,
                   y_title = "population density",
                   x_title = "median income",
                   n_bins = 17) +
  plot_layout(nrow = 1, widths = c(1,.2)) &
  theme(text = element_text(family = "CMU Serif Extra"))