# This code can be used to replicate the figure 1 for "Data alone will not clear the air" paper by E. K. Smith et al. 2025

#### setup #### 
# install.packages("tidyverse")
# install.packages("patchwork")
# install.packages("waffle")
# install.packages("forcats")

library(tidyverse)
library(patchwork)
library(waffle)
library(forcats)

# edit working directory to the folder that contains the dataset
setwd("Smithetal2025")

dir.create("codeOutput")

# load dataset
dataCountry <- read.csv('datasetCommentary.csv', sep = ',', allowEscapes = T)
dataCountry <- subset(dataCountry, !is.na(continent))

##### Panel a) #####
# set up data structure for waffle plot
globalSummary <- dataCountry %>%
  summarize(
    total_referenceGradeCount = sum(referenceGradeCount),
    total_airSensorCount = sum(airSensorCount)
  ) %>% 
  pivot_longer(cols = everything(), 
               names_to = "sensor_type", 
               values_to = "total_count") %>% 
  mutate(sensor_type = factor(sensor_type, 
                              levels = c("total_referenceGradeCount", "total_airSensorCount")))

# Set the number of counts per square
total_squares <- 340  # Fixed total squares
nbrow <- 6

globalSummaryWaffle <- globalSummary

globalSummaryWaffle$total_count <- round((globalSummaryWaffle$total_count / sum(globalSummaryWaffle$total_count)) * total_squares)


# Waffle plot

# tall version for final plot
heights_df_tall <- globalSummary %>%
  mutate(
    cumulative_height = (cumsum(total_count)/180),  # Cumulative total rows for stacking
    midpoint_y = (cumulative_height - ((total_count / 2)/180))/3.3  # Midpoint for each stack
  )

p1_tall <- ggplot(globalSummaryWaffle, aes(fill = sensor_type, values = total_count)) +
  geom_waffle(n_rows = nbrow, size = 1, colour = "white"
              , flip = T
  ) +
  scale_fill_manual(name = NULL,
                    values = c("#BA94CF", "#F6D91B"),
                    labels = c("Regulatory\nGrade", "Non-Regulatory\nGrade")) +
  coord_equal() +
  theme_void()+
  geom_label(
    data = heights_df_tall,
    aes(
      x = 3,  # Fixed x-coordinate
      y = midpoint_y,  # y-coordinate from cumulative calculation
      label = total_count  # Text is raw nb
    ),
    size = 4,  # Font size for the text
    color = "black",  # Text color
    fill = "white",  # Background color of the box
    alpha = 0.7,  # Transparency of the box
    nudge_x = 0.5,
    label.r = unit(0, "pt")
  )


p1 <- ggplot(globalSummaryWaffle, aes(fill = sensor_type, values = total_count)) +
  geom_waffle(n_rows = nbrow, size = 1, colour = "white"
              , flip = F
  ) +
  scale_fill_manual(name = NULL,
                    values = c("#BA94CF", "#F6D91B"),
                    labels = c("Regulatory\nGrade", "Non-Regulatory\nGrade")) +
  coord_equal() +
  theme_void()

p1

# add labels for number of rf and as monitoring stations. finding the right height requires some manual tweaking.
heights_df <- globalSummary %>%
  mutate(
    cumulative_height = (cumsum(total_count)/180),  # Cumulative total rows for stacking
    midpoint_x = (cumulative_height - ((total_count / 2)/180))/3.3  # Midpoint for each stack  # Midpoint for each stack
  )

p1 <- p1 + 
  geom_label(
    data = heights_df,
    aes(
      x = midpoint_x,  # Fixed x-coordinate
      y = 3.5,  # y-coordinate from cumulative calculation
      label = total_count  # Text is raw nb
    ),
    size = 4,  # Font size for the text
    color = "black",  # Text color
    fill = "white",  # Background color of the box
    alpha = 0.7,  # Transparency of the box
    nudge_x = 0.5,
    label.r = unit(0, "pt")
    )

p1

# save the plot
ggsave(plot = p1, filename = paste("codeOutput/panelA.png", sep = ""),
       dpi=600, width = 18, height = 14, units='cm')

##### Panel b) #####

# statistics of each continent
summaryContinent <- dataCountry %>%
  group_by(continent) %>%
  summarize(
    total_referenceGradeCount = sum(referenceGradeCount, na.rm = TRUE),
    total_airSensorCount = sum(airSensorCount, na.rm = TRUE),
    total_allCount = sum(referenceGradeCount, na.rm = TRUE) + sum(airSensorCount, na.rm = TRUE)
  ) %>%
  mutate(
    proportion_referenceGradeGlobal = total_referenceGradeCount / sum(total_referenceGradeCount),
    proportion_airSensorGlobal = total_airSensorCount / sum(total_airSensorCount),
    proportion_referenceGradeIncome = total_referenceGradeCount / (total_referenceGradeCount+total_airSensorCount),
    proportion_airSensorIncome = total_airSensorCount / (total_referenceGradeCount+total_airSensorCount),
  )

# long format for plotting
summaryContinentLong <- summaryContinent %>%
  select(continent, total_referenceGradeCount, total_airSensorCount) %>%
  pivot_longer(cols = contains("total"), names_to = "sensor_type", 
               values_to = "total_count") %>% 
  mutate(sensor_type = factor(sensor_type, 
                              levels = c("total_referenceGradeCount", "total_airSensorCount")))


# quick look at the proportion in EU, NA and Asia compared to the rest
proportionEuNaAsia <- sum(summaryContinent[summaryContinent$continent %in% c("Asia", "North America", "Europe"),]$total_allCount)/sum(summaryContinent$total_allCount)

# reformatting for plotting
summaryContinentLong$continent <- ifelse(
  summaryContinentLong$continent == "North America", "North\nAmerica",
  ifelse(summaryContinentLong$continent == "South America", "South\nAmerica", 
         summaryContinentLong$continent)
)

# change continent order
term_order <- c("Africa",
              "Asia",
              "North\nAmerica",
              "South\nAmerica",
              "Europe",
              "Oceania"
)

summaryContinentLong <- summaryContinentLong %>%
  mutate(continent = factor(continent, levels = term_order))

## Waffle plot p2
summaryContinentLongWaffle <- summaryContinentLong
summaryContinentLongWaffle$total_count <- round((summaryContinentLongWaffle$total_count / sum(summaryContinentLongWaffle$total_count)) * total_squares)

# Create p2
p2 <- ggplot(summaryContinentLongWaffle, aes(fill = sensor_type, values = total_count)) +
  geom_waffle(n_rows = nbrow*2, size = 1, colour = "white"
              , flip = T
  ) +
  scale_fill_manual(name = NULL,
                    values = c("#BA94CF", "#F6D91B"),
                    labels = c("Regulatory\nGrade", "Non-Regulatory\nGrade")) +
  coord_equal() +
  facet_wrap(~continent
             , ncol = 2
             , strip.position = "bottom"
  )+
  theme_void()

p2

# add labels for number of monitoring stations
heights_df <- summaryContinentLong %>%
  group_by(continent) %>% 
  mutate(
    cumulative_height = (cumsum(total_count)/340),  # Cumulative total rows for stacking
    midpoint_y = (cumulative_height - ((total_count / 2)/340))/3.1  # Midpoint for each stack
  )

# need to adjust heights
heights_df <- heights_df %>%
  mutate(
    midpoint_x = case_when(
      continent == "Africa" & sensor_type == "total_referenceGradeCount" ~ 0.5,
      continent == "Africa" & sensor_type != "total_referenceGradeCount" ~ 3.5,
      continent == "South\nAmerica" & sensor_type == "total_referenceGradeCount" ~ 1.5,
      continent == "South\nAmerica" & sensor_type != "total_referenceGradeCount" ~ 5.5,
      continent == "Oceania" & sensor_type == "total_referenceGradeCount" ~ 1,
      continent == "Oceania" & sensor_type != "total_referenceGradeCount" ~ 4.5,
      TRUE ~ 6  # Default value for other continents or conditions
    ),
    midpoint_y = if_else(
      continent %in% c("Africa", "South\nAmerica", "Oceania"),
      2.5,
      midpoint_y
    )
  )

p2 <- p2 + 
  geom_label(
    data = heights_df,
    aes(
      x = midpoint_x,  # Fixed x-coordinate
      y = midpoint_y,  # y-coordinate from cumulative calculation
      label = total_count  # Text is raw nb
    ),
    size = 4,  # Font size for the text
    color = "black",  # Text color
    fill = "white",  # Background color of the box
    alpha = 0.7,  # Transparency of the box
    nudge_x = 0.5,
    label.r = unit(0, "pt")
  )

p2

ggsave(plot = p2, filename = paste("codeOutput/panelB.png", sep = ""),
       dpi=600, width = 18, height = 14, units='cm')


##### Panel c) #####

# statistics of each income group
summaryIncome <- dataCountry %>%
  group_by(Income.group) %>%
  summarize(
    total_referenceGradeCount = sum(referenceGradeCount, na.rm = TRUE),
    total_airSensorCount = sum(airSensorCount, na.rm = TRUE)
  ) %>%
  mutate(
    proportion_referenceGradeGlobal = total_referenceGradeCount / sum(total_referenceGradeCount),
    proportion_airSensorGlobal = total_airSensorCount / sum(total_airSensorCount),
    proportion_referenceGradeIncome = total_referenceGradeCount / (total_referenceGradeCount+total_airSensorCount),
    proportion_airSensorIncome = total_airSensorCount / (total_referenceGradeCount+total_airSensorCount),
  )

# long format for plotting
summaryIncomeLong <- summaryIncome %>%
  select(Income.group, total_referenceGradeCount, total_airSensorCount) %>%
  pivot_longer(cols = contains("total"), names_to = "sensor_type", 
               values_to = "total_count") %>% 
  mutate(sensor_type = factor(sensor_type, 
                              levels = c("total_referenceGradeCount", "total_airSensorCount")))

summaryIncomeLongWaffle <- summaryIncomeLong

summaryIncomeLongWaffle$total_count <- round((summaryIncomeLongWaffle$total_count / sum(summaryIncomeLongWaffle$total_count)) * total_squares)

# reformat for plotting
summaryIncomeLongWaffle$Income.group <- ifelse(
  summaryIncomeLongWaffle$Income.group == "1. High Income", "High\nIncome",
  ifelse(summaryIncomeLongWaffle$Income.group == "2. Upper middle income", "Upper middle\nincome",
         ifelse(summaryIncomeLongWaffle$Income.group == "3. Lower middle income", "Lower middle\nincome",
                ifelse(summaryIncomeLongWaffle$Income.group == "4. Low income", "Low\nincome", summaryIncomeLongWaffle$Income.group)
         )
  )
)
  
# change income level order
term_order_Inc <- c("High\nIncome",
                "Upper middle\nincome",
                "Lower middle\nincome",
                "Low\nincome"
)

summaryIncomeLongWaffle <- summaryIncomeLongWaffle %>%
  mutate(Income.group = factor(Income.group, levels = term_order_Inc))


# create p3

p3 <- ggplot(summaryIncomeLongWaffle, aes(fill = sensor_type, values = total_count)) +
  geom_waffle(n_rows = nbrow*2, size = 1, colour = "white"
              , flip = T
  ) +
  scale_fill_manual(name = NULL,
                    values = c("#BA94CF", "#F6D91B"),
                    labels = c("Regulatory\nGrade", "Non-Regulatory\nGrade")) +
  coord_equal() +
  facet_wrap(~Income.group,
             ncol = 2
             , strip.position = "bottom"
  )+
  theme_void()

p3

# add labels on plot
heights_df <- summaryIncomeLong %>%
  group_by(Income.group) %>% 
  mutate(
    cumulative_height = (cumsum(total_count)/340),  # Cumulative total rows for stacking
    midpoint_y = (cumulative_height - ((total_count / 2)/340))/3  # Midpoint for each stack
  )

# adjust position manually
heights_df <- heights_df %>%
  mutate(
    midpoint_x = case_when(
      Income.group == "3. Lower middle income" & sensor_type == "total_referenceGradeCount" ~ 2,
      Income.group == "3. Lower middle income" & sensor_type != "total_referenceGradeCount" ~ 6.5,
      Income.group == "4. Low income" & sensor_type == "total_referenceGradeCount" ~ 0.3,
      Income.group == "4. Low income" & sensor_type != "total_referenceGradeCount" ~ 3.5,
      TRUE ~ 6  # Default value for other continents or conditions
    ),
    midpoint_y = if_else(
      Income.group %in% c("3. Lower middle income", "4. Low income"),
      2.5,
      midpoint_y
    )
  )

# reformat for plotting
heights_df$Income.group <- ifelse(
  heights_df$Income.group == "1. High Income", "High\nIncome",
  ifelse(heights_df$Income.group == "2. Upper middle income", "Upper middle\nincome",
         ifelse(heights_df$Income.group == "3. Lower middle income", "Lower middle\nincome",
                ifelse(heights_df$Income.group == "4. Low income", "Low\nincome", heights_df$Income.group)
         )
  )
)
  
heights_df <- heights_df %>%
  mutate(Income.group = factor(Income.group, levels = term_order_Inc))

p3 <- p3 + 
  geom_label(
    data = heights_df,
    aes(
      x = midpoint_x,  # Fixed x-coordinate
      y = midpoint_y,  # y-coordinate from cumulative calculation
      label = total_count  # Text is raw nb
    ),
    size = 4,  # Font size for the text
    color = "black",  # Text color
    fill = "white",  # Background color of the box
    alpha = 0.7,  # Transparency of the box
    nudge_x = 0.5,
    label.r = unit(0, "pt")
  )

p3

ggsave(plot = p3, filename = paste("codeOutput/panelC.png", sep = ""),
       dpi=600, width = 18, height = 14, units='cm')


##### Panel d) #####

# statistics of each continent: percentage of people living less than 5km away from a monitor
summaryRatioContinent <- dataCountry %>%
  group_by(continent) %>%
  summarize(ratioRefPop = sum(populationGriddedMonitoredReference, na.rm = TRUE) / sum(populationGriddedTotal, na.rm = TRUE),
            ratioSensorPop = sum(populationGriddedMonitoredAirSensor, na.rm = TRUE) / sum(populationGriddedTotal, na.rm = TRUE),
            ratioAllPop = (sum(populationGriddedMonitoredReference, na.rm = TRUE) + sum(populationGriddedMonitoredAirSensor, na.rm = TRUE)) / sum(populationGriddedTotal, na.rm = TRUE),
            total_referenceGradeCount = sum(referenceGradeCount, na.rm = TRUE),
            total_airSensorCount = sum(airSensorCount, na.rm = TRUE),
            total_allCount = sum(referenceGradeCount, na.rm = TRUE) + sum(airSensorCount, na.rm = TRUE)
            )

# long format for plotting
summaryRatioContinentLong <- summaryRatioContinent %>%
  select(-ratioAllPop) %>%
  pivot_longer(cols = contains("ratio"), names_to = "sensor_type",
               values_to = "ratio") %>%
  mutate(sensor_type = factor(sensor_type,
                              levels = c("ratioRefPop", "ratioSensorPop")))


summaryRatioContinentLong$continent <- ifelse(
  summaryRatioContinentLong$continent == "North America", "North\nAmerica",
  ifelse(summaryRatioContinentLong$continent == "South America", "South\nAmerica",
         summaryRatioContinentLong$continent)
)
# change continent order
term_order <- c("Africa",
                "Asia",
                "North\nAmerica",
                "South\nAmerica",
                "Europe",
                "Oceania"
)

# reorder based on the ratio
summaryRatioContinentLong <- summaryRatioContinentLong %>%
  mutate(continent = fct_reorder(continent, ratio, .fun = sum))

# Calculate cumulative heights
summaryRatioContinentLong <- summaryRatioContinentLong %>%
  group_by(continent) %>%
  mutate(cum_ratio = cumsum(ratio) * 100 - (ratio * 100 / 2)) # Position at the center of each bar segment

summaryRatioContinentLong$cum_ratio <- ifelse(summaryRatioContinentLong$continent == "Africa" & summaryRatioContinentLong$sensor_type == "ratioSensorPop",
                                              summaryRatioContinentLong$cum_ratio + 2, summaryRatioContinentLong$cum_ratio)
# Base ggplot
p4 <- summaryRatioContinentLong %>%
  mutate(sensor_type = factor(sensor_type, levels = c("ratioSensorPop", "ratioRefPop"))) %>%
  ggplot(aes(x = continent, y = ratio * 100, fill = sensor_type)) +
  geom_bar(stat = 'identity', position = 'stack') +
  scale_fill_manual(
    name = NULL,
    values = c("#F6D91B", "#BA94CF"),
    labels = c("Non-Regulatory\nGrade", "Regulatory\nGrade")
  ) +
  labs(
    x = "Continent",
    y = "% Population living <5Km from monitoring station"
  ) +
  ylim(0, 62)+
  theme_minimal()+
  coord_flip()

# Add geom_label for cumulative positioning
p4 <- p4 + geom_label(
  aes(
    x = continent,
    y = cum_ratio,  # Center of each bar segment
    label = sprintf("%.0f%%", ratio * 100)  # Format label as percentage
  ),
  size = 4,  # Font size
  color = "black",  # Text color
  fill = "white",  # Box background
  alpha = 0.7,  # Transparency
  label.r = unit(0, "pt")  # Remove label box corner rounding
)

p4

ggsave(plot = p4, filename = paste("codeOutput/panelD.png", sep = ""),
       dpi=600, width = 20, height = 6, units='cm')

##### Panel e) #####

# statistics of each income level: percentage of people living less than 5km away from a monitor

summaryRatioIncome <- dataCountry %>%
  group_by(Income.group) %>%
  summarize(ratioRefPop = sum(populationGriddedMonitoredReference, na.rm = TRUE) / sum(populationGriddedTotal, na.rm = TRUE),
            ratioSensorPop = sum(populationGriddedMonitoredAirSensor, na.rm = TRUE) / sum(populationGriddedTotal, na.rm = TRUE),
            ratioAllPop = (sum(populationGriddedMonitoredReference, na.rm = TRUE) + sum(populationGriddedMonitoredAirSensor, na.rm = TRUE)) / sum(populationGriddedTotal, na.rm = TRUE),
            total_referenceGradeCount = sum(referenceGradeCount, na.rm = TRUE),
            total_airSensorCount = sum(airSensorCount, na.rm = TRUE),
            total_allCount = sum(referenceGradeCount, na.rm = TRUE) + sum(airSensorCount, na.rm = TRUE)
  )

# long format for plotting
summaryRatioIncomeLong <- summaryRatioIncome %>%
  select(-ratioAllPop) %>%
  pivot_longer(cols = contains("ratio"), names_to = "sensor_type",
               values_to = "ratio") %>%
  mutate(sensor_type = factor(sensor_type,
                              levels = c("ratioRefPop", "ratioSensorPop")))


# reformat for plotting
summaryRatioIncomeLong$Income.group <- ifelse(
  summaryRatioIncomeLong$Income.group == "1. High Income", "High\nIncome",
  ifelse(summaryRatioIncomeLong$Income.group == "2. Upper middle income", "Upper middle\nincome",
         ifelse(summaryRatioIncomeLong$Income.group == "3. Lower middle income", "Lower middle\nincome",
                ifelse(summaryRatioIncomeLong$Income.group == "4. Low income", "Low\nincome", summaryRatioIncomeLong$Income.group)
         )
  )
)
  
# reorder based on ratio
summaryRatioIncomeLong <- summaryRatioIncomeLong %>%
  mutate(Income.group = fct_reorder(Income.group, ratio, .fun = sum))

# Calculate cumulative heights
summaryRatioIncomeLong <- summaryRatioIncomeLong %>%
  group_by(Income.group) %>%
  mutate(cum_ratio = cumsum(ratio) * 100 - (ratio * 100 / 2)) # Position at the center of each segment

summaryRatioIncomeLong$cum_ratio <- ifelse(summaryRatioIncomeLong$Income.group == "Low\nincome" & summaryRatioIncomeLong$sensor_type == "ratioSensorPop",
                                           summaryRatioIncomeLong$cum_ratio + 2, summaryRatioIncomeLong$cum_ratio)

# Base ggplot
p5 <- summaryRatioIncomeLong %>%
  mutate(sensor_type = factor(sensor_type, levels = c("ratioSensorPop", "ratioRefPop"))) %>%
  ggplot(aes(x = Income.group, y = ratio * 100, fill = sensor_type)) +
  geom_bar(stat = 'identity', position = 'stack') +
  scale_fill_manual(
    name = NULL,
    values = c("#F6D91B", "#BA94CF"),
    labels = c("Non-Regulatory\nGrade", "Regulatory\nGrade")
  ) +
  labs(
    x = "Income Group",
    y = "% Population living <5Km from monitoring station"
  ) +
  ylim(0, 62)+
  theme_minimal()+
  coord_flip()

# Add geom_label for cumulative positioning
p5 <- p5 + geom_label(
  aes(
    x = Income.group,
    y = cum_ratio,  # Center of each bar segment
    label = sprintf("%.0f%%", ratio * 100)  # Format label as percentage
  ),
  size = 4,       # Font size
  color = "black", # Text color
  fill = "white",  # Box background color
  alpha = 0.7,     # Transparency of the box
  label.r = unit(0, "pt")  # Remove label box corner rounding
)

p5

ggsave(plot = p5, filename = paste("codeOutput/panelE.png", sep = ""),
       dpi=600, width = 22, height = 6, units='cm')

##### all panels together #####
arranged_plots <- free(p1_tall + ggtitle(label = "a) Global") + theme(
  legend.text = element_text(size = 14),
  plot.title = element_text(face = 'bold', size = 16),
  legend.key.height = unit(1, "cm"),
  legend.key.width = unit(1, "cm"),
  legend.position="right"
)) | 
  free(p2 + ggtitle(label = "b) Split by Region\n") + theme(
    legend.position = "none",
    strip.text = element_text(size = 16),
    plot.title = element_text(face = 'bold', size = 16)
  )) |
  free(p3 + ggtitle(label = "c) Split by Income\n\n\n") + theme(
    legend.position = "none",
    strip.text = element_text(size = 16),
    plot.title = element_text(face = 'bold', size = 16)
  )) |
  free(p4 + ggtitle(label = "d) Split by Region\n") + theme(
    legend.position = "none",
    axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_blank(),
    # axis.title.y = element_blank(),
    axis.title.y = element_text(size = 14),
    strip.text = element_text(size = 16),
    plot.title = element_text(face = 'bold', size = 16)
  )) |
  free(p5 + ggtitle(label = "e) Split by Income\n") + theme(
    legend.position = "none",
    axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    # axis.title.y = element_text(size = 14),
    strip.text = element_text(size = 16),
    plot.title = element_text(face = 'bold', size = 16)
  ))

ggsave("codeOutput/Figure1_waffle_horizontal_2.pdf", arranged_plots, width = 19, height = 9)
ggsave("codeOutput/Figure1_waffle_horizontal_2.png", arranged_plots, width = 19, height = 9)