Trumpet plot

gistfile1.txt

library(tidyverse)
library(ggsankey)

## data via https://ourworldindata.org/grapher/total-agricultural-area-over-the-long-term
hablar::set_wd_to_script_path()


df <-
  read_csv("total-agricultural-area-over-the-long-term.csv") %>%
  janitor::clean_names() %>%
  filter(!entity %in% c("World", "Greenland"), year >= 0) %>%
  mutate(entity = case_when(
    #str_detect(entity, "Brazil|United States|Canada") ~ "Americas and the Carribean",
    str_detect(entity, "Brazil") ~ "Latin America and\nthe Caribbean",
    str_detect(entity, "Canada|United States") ~ "Northern America",
    str_detect(entity, "Russia") ~ "Europe",
    str_detect(entity, "Rest of Asia") ~ "Asia (excluding\nIndia and China)",
    #str_detect(entity, "India|China") ~ "Asia",
    TRUE ~ entity
  ))

df <- df %>%
  group_by(entity, year) %>%
  summarize(cropland = sum(agricultural_area_crops_grazing_hyde_2017, na.rm = TRUE)) %>%
  #mutate(timestep = year) %>%
  #mutate(timestep = (year %/% 100) * 100) %>%
  mutate(timestep = ((year * 2) %/% 100) * 100 / 2) %>%
  group_by(entity, timestep) %>%
  summarize(cropland = max(cropland, na.rm = TRUE)) %>%
  group_by(entity) %>%
  mutate(last = cropland[which(timestep == max(timestep))]) %>%
  ungroup() %>%
  mutate(entity = fct_reorder(entity, -last))

df_labs <-
  tibble(
    entity = factor(levels(df$entity), levels = levels(df$entity)),
    cropland =  c(720000000, 400000000, 125000000, -125000000, -340000000,
                  -540000000, -710000000, -815000000, -872000000)
  )

cols <- c("firebrick", "tan", "#F39530", "#778e99", "grey74",
          "#003b4c", "#bbcad2", "grey81", "grey88")

set.seed(1)

###
row_shift_arc <- function(xmiddle, direction, ymin1, ymax1, ymin2, ymax2, radius, thickness) {

  # inner arc
  thickness_mid <- ((ymax1 - ymin1) + (ymax2 - ymin2))/2

  data <- tibble(x0 = xmiddle, y0 = ymax1 + (ymin2 - ymax1)/2, a = radius, b = (ymin2 - ymax1)/2, angle = 0)
  data$m1 <- 2
  data$m2 <- 2
  n_ellipses <- nrow(data)
  n <- 10000

  data <- data[rep(seq_len(n_ellipses), each = n), ]
  points <- rep(seq(0, 2 * pi, length.out = n + 1)[seq_len(n)],
                n_ellipses)
  cos_p <- cos(points)
  sin_p <- sin(points)
  x_tmp <- abs(cos_p)^(2 / data$m1) * data$a * sign(cos_p)
  y_tmp <- abs(sin_p)^(2 / data$m2) * data$b * sign(sin_p)
  data$x <- data$x0 + x_tmp * cos(data$angle) - y_tmp * sin(data$angle)
  data$y <- data$y0 + x_tmp * sin(data$angle) + y_tmp * cos(data$angle)

  inner_arc <- data

  # outer arc
  data <- tibble(x0 = xmiddle, y0 = ymin1 + (ymax2 - ymin1)/2, a = radius + thickness, b = (ymax2 - ymin1)/2, angle = 0)
  data$m1 <- 2
  data$m2 <- 2
  n_ellipses <- nrow(data)
  n <- 10000

  data <- data[rep(seq_len(n_ellipses), each = n), ]
  points <- rep(seq(0, 2 * pi, length.out = n + 1)[seq_len(n)],
                n_ellipses)
  cos_p <- cos(points)
  sin_p <- sin(points)
  x_tmp <- abs(cos_p)^(2 / data$m1) * data$a * sign(cos_p)
  y_tmp <- abs(sin_p)^(2 / data$m2) * data$b * sign(sin_p)
  data$x <- data$x0 + x_tmp * cos(data$angle) - y_tmp * sin(data$angle)
  data$y <- data$y0 + x_tmp * sin(data$angle) + y_tmp * cos(data$angle)

  outer_arc <- data

  data <- bind_rows(
    inner_arc %>% arrange(y),
    outer_arc %>% arrange(-y)
  )

  if(direction == "right") {
    data <- data %>%
      filter(x >= xmiddle)
  }

  if(direction == "left") {
    data <- data %>%
      filter(x <= xmiddle)
  }

  data
}

row_shift_arc(500, "right", -10, 0, 150, 160, 30) %>%
  ggplot(aes(x, y)) +
  geom_polygon() +
  coord_fixed()



# Test
space_var = 30e6

df1 <- df %>% filter(timestep <= 500) %>% mutate(shift = 0e08, timestep = timestep - 0)
df2 <- df %>% filter(timestep > 500, timestep <= 1000) %>% mutate(shift = 10e08, timestep = timestep - 500, timestep = abs(timestep - 600))
df3 <- df %>% filter(timestep > 1000) %>% mutate(shift = 2 * 10e08, timestep = timestep - 1000)


# Arc 1
tt <- df1 %>%
  filter(timestep == max(timestep)) %>%
  arrange(cropland) %>%
  mutate(total = sum(cropland) + (n() - 1) * space_var,
         total_max = shift - total/2,
         ymax1 = total_max + cumsum(cropland) + (row_number() - 1) * space_var,
         ymin1 = ymax1 - cropland) %>%
  select(entity, timestep, ymax1, ymin1)

tt2 <- df2 %>%
  filter(timestep == max(timestep)) %>%
  arrange(cropland) %>%
  mutate(total = sum(cropland) + (n() - 1) * space_var,
         total_max = shift - total/2,
         ymax2 = total_max + cumsum(cropland) + (row_number() - 1) * space_var,
         ymin2 = ymax2 - cropland) %>%
  select(entity, timestep, ymax2, ymin2)

arcs <- tt %>%
  left_join(tt2, by = c("entity", "timestep"))

arc_data <- map_dfr(arcs$entity %>% unique(),
        ~{
          df <- arcs %>%
            filter(entity == .x)

          row_shift_arc(df$timestep, "right", df$ymin1, df$ymax1, df$ymin2, df$ymax2, 100, 20) %>%
            mutate(entity = .x)
        })

# Arc 2
tt <- df2 %>%
  filter(timestep == min(timestep)) %>%
  arrange(cropland) %>%
  mutate(total = sum(cropland) + (n() - 1) * space_var,
         total_max = shift - total/2,
         ymax1 = total_max + cumsum(cropland) + (row_number() - 1) * space_var,
         ymin1 = ymax1 - cropland) %>%
  select(entity, timestep, ymax1, ymin1)

tt2 <- df3 %>%
  filter(timestep == min(timestep)) %>%
  arrange(cropland) %>%
  mutate(total = sum(cropland) + (n() - 1) * space_var,
         total_max = shift - total/2,
         ymax2 = total_max + cumsum(cropland) + (row_number() - 1) * space_var,
         ymin2 = ymax2 - cropland) %>%
  select(entity, timestep, ymax2, ymin2)

arcs <- tt %>%
  left_join(tt2, by = c("entity", "timestep"))

arc_data2 <- map_dfr(arcs$entity %>% unique(),
                    ~{
                      df <- arcs %>%
                        filter(entity == .x)

                      row_shift_arc(df$timestep, "left", df$ymin1, df$ymax1, df$ymin2, df$ymax2, 100, 20) %>%
                        mutate(entity = .x)
                    })

# Plot
ggplot(bind_rows(df1, df2),
       aes(x = timestep, value = cropland, node = entity, fill = entity, color = entity)) +
  geom_sankey_bump(data = df1, aes(shift = shift), space = space_var, smooth = 6) +
  geom_sankey_bump(data = df2, aes(shift = shift), space = space_var, smooth = 6) +
  geom_sankey_bump(data = df3, aes(shift = shift), space = space_var, smooth = 6) +
  geom_polygon(data = arc_data, aes(x, y, fill = entity, group = entity, color = entity), inherit.aes = FALSE, alpha = 1) +
  geom_polygon(data = arc_data2, aes(x, y, fill = entity, group = entity, color = entity), inherit.aes = FALSE, alpha = 1) +
  # geom_text(
  #   data = df_labs,
  #   aes(x = 2003, y = cropland, label = entity,
  #       color = entity, color = after_scale(colorspace::darken(color, .2))),
  #   hjust = 0, size = 6.5, family = "Bebas Neue", lineheight = .75
  # ) +
  # scale_x_continuous(expand = c(0, 0), limits = c(NA, 2100),
  #                    breaks = c(500, 1000, 1500, 1600, 1700, 1800, 1900)) +
  # scale_y_continuous(expand = c(.01, .01)) +
  scale_color_manual(values = cols, guide = FALSE) +
  scale_fill_manual(values = cols, guide = FALSE) +
  theme_void() +
theme(plot.background = element_rect(fill = "#FAF5F0"),
      panel.grid.major.x = element_line(color = "#E8E1DB")) +
NULL