# Séance 5 : Visualiser avec ggplot2 (et cartes)
# Source : https://pierrebeaucoral.github.io/r-for-beginners/fr/
# Ouvrez votre projet r-course dans RStudio, puis exécutez ce fichier ligne par ligne.

Sys.setenv(LANGUAGE = "en")   # messages de R en anglais

library(tidyverse)
library(scales)
library(ggrepel)
library(sf)
library(rnaturalearth)
library(here)

source(here("R", "theme_cours.R"))   # theme_cours(), pal_cours, pal_regions

wdi <- read_csv(here("data", "wdi_panel.csv"), show_col_types = FALSE)

# Une ligne par pays pour une année : la coupe transversale utilisée plus bas
wdi_2022 <- wdi |>
  filter(year == 2022)

wdi_2022 |>
  drop_na(life_exp) |>
  slice_min(life_exp, n = 3) |>
  select(country, iso3c, life_exp)

# L'espérance de vie de la RCA n'est pas fiable dans cette version des WDI (14,7 en 2009, 18,8 en 2022) : on l'exclut pour toutes les années
wdi <- wdi |> mutate(life_exp = if_else(iso3c == "CAF", NA_real_, life_exp))

wdi_2022 <- filter(wdi, year == 2022)

wdi_2022_ok <- wdi_2022 |>
  drop_na(gdp_pc_ppp, life_exp)        # un point a besoin des deux coordonnées

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp))

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point()

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point() +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale()))

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
  geom_point(alpha = 0.8) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  scale_colour_manual(values = pal_regions, guide = "none")

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
  geom_point(alpha = 0.8) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  scale_colour_manual(values = pal_regions) +
  guides(colour = guide_legend(nrow = 4, title.position = "top")) +
  labs(
    title = str_wrap("Les pays plus riches vivent plus longtemps, mais les gains s'aplatissent en haut", 60),
    x = "PIB par habitant, PPA (en $ internationaux constants de 2021, échelle log)",
    y = "Espérance de vie à la naissance (ans)",
    colour = "Région",
    caption = "Source : Banque mondiale, World Development Indicators (WDI), 2022."
  ) +
  theme_cours(base_size = 11)

focus <- c("China", "India", "Nigeria", "United States")

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point(aes(size = population, colour = region), alpha = 0.7) +
  geom_text_repel(
    data = filter(wdi_2022_ok, country %in% focus),
    aes(label = country),
    min.segment.length = 0, box.padding = 1, point.padding = 0.5, seed = 1, size = 4.2,
    bg.color = "white", bg.r = 0.15
  ) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  scale_size_area(max_size = 14, breaks = c(1e7, 1e8, 1e9),
                  labels = c("10 M", "100 M", "1 B")) +
  scale_colour_manual(values = pal_regions) +
  guides(colour = guide_legend(nrow = 3, title.position = "top", order = 1,
                               override.aes = list(size = 4)),
         size = guide_legend(title.position = "top", order = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical",
        legend.box.just = "left", legend.spacing.y = unit(2, "pt")) +
  labs(
    title = str_wrap("Revenu et espérance de vie progressent ensemble, mais des revenus proches peuvent cacher 20 ans d'écart", 60),
    x = "PIB par habitant, PPA (en $ internationaux constants de 2021, échelle log)",
    y = "Espérance de vie à la naissance (ans)",
    colour = "Région", size = "Population",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Un point par pays ; les pays aux données manquantes sont exclus.", 95)
  )

wdi_2022 |>
  drop_na(elec_access, life_exp, income) |>
  ggplot(aes(x = elec_access, y = life_exp, colour = income)) +
  geom_point(alpha = 0.8, size = 2) +
  scale_x_continuous(labels = label_number(suffix = "%")) +
  scale_colour_manual(values = c("Low income" = "#D55E00", "Lower middle income" = "#E69F00",
                                 "Upper middle income" = "#56B4E9", "High income" = "#0072B2")) +
  guides(colour = guide_legend(nrow = 2)) +
  labs(
    title = str_wrap("Les pays sans électricité universelle ont une espérance de vie nettement plus basse", 60),
    x = "Population ayant accès à l'électricité (%)",
    y = "Espérance de vie à la naissance (ans)",
    colour = str_wrap("Groupe de revenu (classification actuelle de la Banque mondiale)", 25),
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022.", 95)
  )

compare <- c("India", "China", "Nigeria", "Brazil", "United States", "Germany", "Japan", "Bangladesh")
highlight <- "India"

life_lines <- wdi |>
  filter(country %in% compare) |>
  drop_na(life_exp)

ggplot(life_lines, aes(x = year, y = life_exp, group = country)) +
  geom_line(data = filter(life_lines, country != highlight),
            colour = pal_cours[["gris"]], linewidth = 0.8) +
  geom_line(data = filter(life_lines, country == highlight),
            colour = pal_cours[["accent"]], linewidth = 1.6) +
  geom_text_repel(data = filter(life_lines, year == max(year)),
                  aes(label = country, colour = country == highlight),
                  hjust = 0, direction = "y", xlim = c(2023.4, NA), size = 4,
                  segment.color = NA, show.legend = FALSE, seed = 1) +
  scale_colour_manual(values = c(`TRUE` = pal_cours[["accent"]], `FALSE` = "grey40")) +
  scale_x_continuous(breaks = seq(2000, 2020, 5), expand = expansion(mult = c(0.02, 0.02))) +
  coord_cartesian(xlim = c(2000, 2027), clip = "off") +
  labs(
    title = str_wrap("L'Inde a gagné neuf ans d'espérance de vie depuis 2000 ; parmi ces huit pays, seul le Bangladesh a gagné davantage", 60),
    x = "Année",
    y = "Espérance de vie à la naissance (ans)",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2000-2023. Inde en rouge ; sept pays de comparaison en gris.", 95)
  )

life_region <- wdi |>
  drop_na(life_exp) |>
  summarise(life_med = median(life_exp), .by = c(region, year))

ggplot(life_region, aes(x = year, y = life_med, colour = region)) +
  geom_line(linewidth = 1.1) +
  geom_text_repel(
    data = filter(life_region, year == max(year)),
    aes(label = str_wrap(region, 20)),
    hjust = 0, direction = "y", xlim = c(2023.4, NA), size = 3.4, lineheight = 0.85,
    segment.color = NA, show.legend = FALSE
  ) +
  scale_colour_manual(values = pal_regions, guide = "none") +
  scale_x_continuous(breaks = seq(2000, 2020, 5), expand = expansion(mult = c(0.02, 0.05))) +
  coord_cartesian(xlim = c(2000, 2028), clip = "off") +
  labs(
    title = str_wrap("L'Afrique subsaharienne a gagné 12 ans depuis 2000 mais reste derrière toutes les autres régions", 60),
    x = "Année",
    y = "Espérance de vie médiane à la naissance (ans)",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2000-2023. Médiane des pays de chaque région (pays non pondérés par la population).", 95)
  )

gdp_region <- wdi_2022_ok |>
  summarise(gdp_med = median(gdp_pc_ppp), n_countries = n(), .by = region)

gdp_region |>
  mutate(region = fct_reorder(region, gdp_med)) |>
  ggplot(aes(x = gdp_med, y = region)) +
  geom_col(fill = pal_cours[["principal"]], width = 0.7) +
  geom_text(aes(label = label_dollar(accuracy = 0.1, scale_cut = cut_short_scale())(gdp_med)), hjust = -0.1, size = 4) +
  scale_x_continuous(labels = label_dollar(scale_cut = cut_short_scale()), expand = expansion(mult = c(0, 0.15))) +
  scale_y_discrete(labels = label_wrap(28)) +
  labs(
    title = str_wrap("Un pays nord-américain typique est environ sept fois plus riche qu'un pays d'Asie du Sud", 60),
    x = "PIB par habitant médian, PPA (en $ internationaux constants de 2021)",
    y = NULL,
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Médiane des pays ; l'Amérique du Nord ne compte que 3 pays.", 95)
  )

niveaux <- c("Low income", "Lower middle income", "Upper middle income", "High income")

wdi_2022 |>
  mutate(income = factor(income, levels = niveaux, ordered = TRUE)) |>
  ggplot(aes(x = income)) +
  geom_bar(fill = pal_cours[["principal"]], width = 0.7) +
  geom_text(aes(label = after_stat(count)), stat = "count", vjust = -0.4, size = 4.5) +
  scale_x_discrete(labels = label_wrap(14)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  labs(
    title = str_wrap("Les pays à revenu élevé forment le plus grand groupe, mais 25 pays restent à faible revenu", 60),
    x = "Groupe de revenu (classification actuelle de la Banque mondiale)",
    y = "Nombre de pays",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Toutes les économies du panel WDI.", 95)
  )

wdi_2022_ok |>
  slice_max(gdp_pc_ppp, n = 10) |>
  mutate(country = fct_reorder(country, gdp_pc_ppp)) |>
  ggplot(aes(x = gdp_pc_ppp, y = country, fill = region)) +
  geom_col(width = 0.7) +
  scale_x_continuous(labels = label_dollar(scale_cut = cut_short_scale()), expand = expansion(mult = c(0, 0.05))) +
  scale_fill_manual(values = pal_regions) +
  guides(fill = guide_legend(nrow = 4)) +
  labs(
    title = str_wrap("Neuf des dix pays les plus riches sont de petites économies ou des places financières ; les États-Unis font exception", 60),
    x = "PIB par habitant, PPA (en $ internationaux constants de 2021)",
    y = NULL, fill = "Région",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022.", 95)
  )

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp)) +
  geom_histogram(aes(y = after_stat(density)), bins = 25,
                 fill = pal_cours[["principal"]], alpha = 0.6, colour = "white") +
  geom_density(colour = pal_cours[["accent"]], linewidth = 1) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = str_wrap("Les revenus mondiaux s'étalent sur plus d'un facteur 100, sans pays typique unique", 60),
    x = "PIB par habitant, PPA (en $ internationaux constants de 2021, échelle log)",
    y = "Densité",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Une observation par pays ; courbe rouge = estimation de densité par noyau.", 95)
  )

wdi_2022_ok |>
  mutate(income = factor(income, levels = niveaux, ordered = TRUE)) |>
  ggplot(aes(x = income, y = gdp_pc_ppp)) +
  geom_boxplot(outliers = FALSE, fill = "grey92", colour = "grey30", width = 0.55) +
  geom_jitter(width = 0.15, alpha = 0.5, colour = pal_cours[["principal"]], size = 1.6) +
  scale_x_discrete(labels = label_wrap(14)) +
  scale_y_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = str_wrap("Le revenu médian augmente fortement d'un groupe de revenu au suivant, mais les groupes voisins se chevauchent", 60),
    x = "Groupe de revenu (classification actuelle de la Banque mondiale)",
    y = str_wrap("PIB par habitant, PPA (en $ internationaux constants de 2021, échelle log)", 40),
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Boîte = quartiles et médiane ; points = pays.", 95)
  )

wdi |>
  filter(year %in% c(2000, 2022)) |>
  drop_na(life_exp) |>
  ggplot(aes(x = life_exp)) +
  geom_histogram(binwidth = 2.5, fill = pal_cours[["principal"]], colour = "white") +
  facet_wrap(~ year) +
  labs(
    title = str_wrap("La queue des faibles espérances de vie s'est réduite entre 2000 et 2022", 60),
    x = "Espérance de vie à la naissance (ans)",
    y = "Nombre de pays",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2000 et 2022. Une observation par pays et par année.", 95)
  )

world <- ne_countries(scale = "medium", returnclass = "sf") |>
  filter(name != "Antarctica")          # pas de données là-bas, et il mange un cinquième de la carte

dim(world)
world |>
  st_drop_geometry() |>
  select(name, iso_a3, iso_a3_eh, adm0_a3) |>
  slice(1:5)

# Pays WDI sans correspondance dans la carte, pour chaque code candidat
unmatched <- function(code) {
  wdi_2022 |>
    anti_join(st_drop_geometry(world), by = c("iso3c" = code)) |>
    pull(country)
}

# Combien d'économies WDI restent sans correspondance avec chaque code ?
length(unmatched("iso_a3"))
length(unmatched("iso_a3_eh"))
length(unmatched("adm0_a3"))

unmatched("iso_a3_eh")

world |>
  st_drop_geometry() |>
  anti_join(wdi_2022, by = c("iso_a3_eh" = "iso3c")) |>
  select(name, iso_a3_eh) |>
  head(8)

world_gdp <- world |>
  select(iso_a3_eh, geometry) |>
  left_join(select(wdi_2022, iso3c, country, gdp_pc_ppp, elec_access),
            by = c("iso_a3_eh" = "iso3c"))

sum(is.na(world_gdp$gdp_pc_ppp))   # formes dessinées en gris (pas de données)

cuts   <- c(0, 2500, 5000, 10000, 20000, 40000, Inf)
labs_g <- c("< $2,500", "$2,500-5,000", "$5,000-10,000", "$10,000-20,000",
            "$20,000-40,000", "> $40,000")

world_gdp |>
  mutate(gdp_class = cut(gdp_pc_ppp, breaks = cuts, labels = labs_g, right = FALSE)) |>
  ggplot() +
  geom_sf(aes(fill = fct_na_value_to_level(gdp_class, "Pas de données")),   # NA devient un niveau
          colour = "white", linewidth = 0.1) +
  scale_fill_manual(values = c(hcl.colors(6, "YlGnBu", rev = TRUE), "grey85"),  # une couleur par niveau, dans l'ordre
                    drop = FALSE) +
  coord_sf(crs = "+proj=robin") +
  guides(fill = guide_legend(nrow = 2, title.position = "top")) +
  labs(
    x = NULL, y = NULL,
    title = str_wrap("Le revenu par habitant est le plus élevé au Nord et le plus bas en Afrique subsaharienne", 60),
    fill = "PIB par habitant, PPA (en $ internationaux constants de 2021)",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Gris clair : pas de données.", 95)
  ) +
  theme(axis.text = element_blank(), panel.grid.major = element_line(colour = "grey94"))

bottom3 <- world_gdp |>
  drop_na(elec_access) |>
  slice_min(elec_access, n = 3, with_ties = FALSE) |>
  mutate(label = paste0(str_remove(country, ",.*"), " ", round(elec_access), "%")) |>
  st_point_on_surface()                       # un point à l'intérieur de chaque polygone

bottom3$label                                  # vérifier quels pays sont annotés

world_gdp |>
  mutate(elec_class = cut(elec_access, breaks = c(0, 50, 80, 99, 100.01),
                          labels = c("< 50%", "50-80%", "80-99%", "99-100%"),
                          right = FALSE)) |>
  ggplot() +
  geom_sf(aes(fill = fct_na_value_to_level(elec_class, "Pas de données")),
          colour = "white", linewidth = 0.1) +
  geom_label_repel(data = bottom3, aes(label = label, geometry = geometry),
                   stat = "sf_coordinates", size = 3.4, fill = alpha("white", 0.85),
                   min.segment.length = Inf, box.padding = 0.2, seed = 1) +  # pas de traits de repère
  scale_fill_manual(values = c("#08306b", "#2171b5", "#9ecae1", "#eff3ff", "grey85"),  # dans l'ordre des niveaux
                    drop = FALSE) +
  coord_sf(crs = "+proj=robin") +
  guides(fill = guide_legend(nrow = 2, title.position = "top")) +
  labs(
    x = NULL, y = NULL,
    title = str_wrap("20 des 22 pays où moins de la moitié de la population a l'électricité sont en Afrique subsaharienne", 60),
    fill = "Population ayant accès à l'électricité",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Étiquettes (placées près des pays) : les trois pays au plus faible accès. Gris clair : pas de données.", 95)
  ) +
  theme(axis.text = element_blank(), panel.grid.major = element_line(colour = "grey94"))

# dir.create(here("output"), showWarnings = FALSE)   # output/ est créé s'il manque
# 
# p_gapminder <- ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
#   geom_point() +
#   scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale()))
# 
# ggsave(here("output", "gapminder.png"), plot = p_gapminder,
#        width = 8, height = 5, units = "in", dpi = 300)    # pour Word / diapositives
# ggsave(here("output", "gapminder.pdf"), plot = p_gapminder,
#        width = 8, height = 5)                             # vectoriel, pour LaTeX

p_life_region <- wdi_2022 |>
  drop_na(life_exp) |>
  mutate(region = fct_reorder(region, life_exp, .fun = median)) |>
  ggplot(aes(x = life_exp, y = region)) +
  geom_boxplot(fill = "grey92", colour = "grey30", width = 0.55) +
  scale_y_discrete(labels = label_wrap(28)) +
  labs(
    title = str_wrap("L'Afrique subsaharienne a l'espérance de vie médiane la plus basse, environ 9 ans sous la région suivante", 60),
    x = "Espérance de vie à la naissance (ans)",
    y = NULL,
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2022. Une observation par pays.", 95)
  )

p_life_region

# dir.create(here("output"), showWarnings = FALSE)
# ggsave(here("output", "life_exp_region.png"), plot = p_life_region,
#        width = 8, height = 5, dpi = 300)