Mozambique January 2026 Flood Event

library(tidyverse)
library(cumulus)
library(DBI)
library(zoo)
library(dbplyr)

con <- pg_con()
df_lookup <- cumulus::blob_load_admin_lookup()

provinces <- c("Maputo", "Sofala", "Gaza")

df_lookup_aoi <- df_lookup |>
  filter(
    ADM_LEVEL == 1,
    ISO3 == "MOZ",
    ADM1_NAME %in% provinces
  )

imerg <- tbl(con, "imerg") |>
  filter(
    adm_level == 1,
    pcode %in% df_lookup_aoi$ADM1_PCODE
  ) |>
  collect() |>
  left_join(
    df_lookup_aoi |> select(pcode = ADM1_PCODE, adm1_name = ADM1_NAME),
    by = "pcode"
  ) |>
  select(iso3, pcode, adm1_name, valid_date, mean) |>
  arrange(adm1_name, valid_date)

Rainfall Analysis (IMERG)

Recent Rainfall Context

Method: Daily IMERG satellite precipitation aggregated to province level.

Key finding: Mid-January 2026 shows a distinct rainfall spike across all three provinces, with peak daily totals reaching 30-50+ mm/day.

imerg_recent <- imerg |>
  filter(valid_date >= max(valid_date) - 60)

top_3_days <- imerg_recent |>
  group_by(adm1_name) |>
  slice_max(mean, n = 3) |>
  ungroup()

ggplot(imerg_recent, aes(x = valid_date, y = mean, color = adm1_name)) +
  geom_line(linewidth = 0.6) +
  geom_point(size = 1) +
  geom_point(data = top_3_days, size = 2.5, shape = 21, fill = "white", stroke = 1.2) +
  geom_text(
    data = top_3_days,
    aes(label = format(valid_date, "%b %d")),
    vjust = -1.2, size = 2.5, show.legend = FALSE
  ) +
  facet_wrap(~adm1_name, ncol = 1) +
  scale_x_date(date_breaks = "2 weeks", date_labels = "%b %d") +
  labs(
    title = "Daily Precipitation - Last 60 Days",
    subtitle = "Top 3 days labeled per province",
    x = NULL,
    y = "Precipitation (mm/day)",
    caption = "Data: IMERG"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    strip.text = element_text(face = "bold", size = 11),
    panel.grid.minor = element_blank()
  )

January 2026 Accumulations

Method: Rolling 3-day and 7-day precipitation sums identify peak accumulation windows.

Key finding: Maximum accumulations occurred mid-January. Shaded regions show the worst 3-day and 7-day periods per province.

imerg_rolling <- imerg |>
  group_by(adm1_name) |>
  mutate(
    roll3_sum = rollsum(mean, k = 3, fill = NA, align = "right"),
    roll7_sum = rollsum(mean, k = 7, fill = NA, align = "right")
  ) |>
  ungroup()

jan_2026 <- imerg_rolling |>
  filter(year(valid_date) == 2026, month(valid_date) == 1)

max_3day <- jan_2026 |>
  filter(!is.na(roll3_sum)) |>
  group_by(adm1_name) |>
  slice_max(roll3_sum, n = 1) |>
  mutate(
    start_date = valid_date - 2,
    period_label = paste0(format(start_date, "%b %d"), "-", format(valid_date, "%d"))
  ) |>
  ungroup()

max_7day <- jan_2026 |>
  filter(!is.na(roll7_sum)) |>
  group_by(adm1_name) |>
  slice_max(roll7_sum, n = 1) |>
  mutate(
    start_date = valid_date - 6,
    period_label = paste0(format(start_date, "%b %d"), "-", format(valid_date, "%d"))
  ) |>
  ungroup()

ggplot(jan_2026, aes(x = valid_date, y = mean)) +
  geom_rect(
    data = max_7day,
    aes(xmin = start_date, xmax = valid_date, ymin = -Inf, ymax = Inf),
    fill = "steelblue", alpha = 0.15, inherit.aes = FALSE
  ) +
  geom_rect(
    data = max_3day,
    aes(xmin = start_date, xmax = valid_date, ymin = -Inf, ymax = Inf),
    fill = "steelblue", alpha = 0.25, inherit.aes = FALSE
  ) +
  geom_col(fill = "steelblue", width = 0.8) +
  geom_label(
    data = max_3day,
    aes(x = start_date + 1, y = Inf,
        label = paste0("3-day: ", round(roll3_sum, 0), " mm")),
    vjust = 1.2, hjust = 0, size = 3.5, fontface = "bold",
    color = "grey20", fill = "white", label.size = 0.3, label.r = unit(0.15, "lines")
  ) +
  geom_label(
    data = max_7day,
    aes(x = start_date + 3, y = Inf,
        label = paste0("7-day: ", round(roll7_sum, 0), " mm")),
    vjust = 2.8, hjust = 0, size = 3.5,
    color = "grey40", fill = "white", label.size = 0.3, label.r = unit(0.15, "lines")
  ) +
  facet_wrap(~adm1_name, ncol = 1, scales = "free_y") +
  scale_x_date(date_breaks = "1 week", date_labels = "%b %d") +
  labs(
    title = "January 2026 Daily Precipitation",
    subtitle = "Shaded: max 3-day (dark) and 7-day (light) accumulation periods",
    x = NULL,
    y = "Precipitation (mm/day)"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold", size = 11),
    panel.grid.minor = element_blank()
  )

Return Period Analysis

Method: Empirical return periods using Weibull plotting position. January 2026 event compared against historical annual maxima from IMERG record.

Key finding: The 7-day accumulation shows higher return periods than the 3-day, indicating sustained rainfall over a full week was the more exceptional aspect of this event.

3-Day Sum

annual_max_3day <- imerg_rolling |>
  filter(!is.na(roll3_sum)) |>
  mutate(year = year(valid_date)) |>
  group_by(adm1_name, year) |>
  slice_max(roll3_sum, n = 1) |>
  ungroup() |>
  group_by(adm1_name) |>
  mutate(
    n_years = n(),
    rank = rank(-roll3_sum, ties.method = "first"),
    return_period = (n_years + 1) / rank
  ) |>
  ungroup()

rp_3day <- max_3day |>
  left_join(
    annual_max_3day |> select(adm1_name, annual_max = roll3_sum),
    by = "adm1_name", relationship = "many-to-many"
  ) |>
  group_by(adm1_name, start_date, valid_date, roll3_sum, period_label) |>
  summarise(
    n_exceeded = sum(roll3_sum > annual_max),
    n_years = n(),
    rp_years = (n_years + 1) / pmax(1, n_years - n_exceeded),
    .groups = "drop"
  )

top10_3day <- annual_max_3day |>
  group_by(adm1_name) |>
  slice_min(rank, n = 10) |>
  ungroup()

ggplot(annual_max_3day, aes(x = return_period, y = roll3_sum)) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 3, alpha = 0.5) +
  geom_text(
    data = top10_3day,
    aes(label = paste0("'", substr(year, 3, 4))),
    hjust = -0.5, vjust = 0.5, size = 2.5
  ) +
  geom_hline(
    data = rp_3day,
    aes(yintercept = roll3_sum),
    color = "#D62828", linetype = "dashed", linewidth = 1
  ) +
  geom_label(
    data = rp_3day,
    aes(x = Inf, y = roll3_sum,
        label = paste0(period_label, "\n~", round(rp_years, 1), "-yr")),
    hjust = 1.05, vjust = 0.5, size = 2.5,
    label.size = 0, fill = "white", color = "#D62828", fontface = "bold"
  ) +
  scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50)) +
  facet_wrap(~adm1_name, scales = "free_y") +
  labs(
    title = "Return Period: 3-Day Sum",
    subtitle = "January 2026 event in red",
    x = "Return Period (years)",
    y = "3-Day Sum (mm)"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold", size = 11),
    panel.grid.minor = element_blank()
  )

7-Day Sum

annual_max_7day <- imerg_rolling |>
  filter(!is.na(roll7_sum)) |>
  mutate(year = year(valid_date)) |>
  group_by(adm1_name, year) |>
  slice_max(roll7_sum, n = 1) |>
  ungroup() |>
  group_by(adm1_name) |>
  mutate(
    n_years = n(),
    rank = rank(-roll7_sum, ties.method = "first"),
    return_period = (n_years + 1) / rank
  ) |>
  ungroup()

rp_7day <- max_7day |>
  left_join(
    annual_max_7day |> select(adm1_name, annual_max = roll7_sum),
    by = "adm1_name", relationship = "many-to-many"
  ) |>
  group_by(adm1_name, start_date, valid_date, roll7_sum, period_label) |>
  summarise(
    n_exceeded = sum(roll7_sum > annual_max),
    n_years = n(),
    rp_years = (n_years + 1) / pmax(1, n_years - n_exceeded),
    .groups = "drop"
  )

top10_7day <- annual_max_7day |>
  group_by(adm1_name) |>
  slice_min(rank, n = 10) |>
  ungroup()

ggplot(annual_max_7day, aes(x = return_period, y = roll7_sum)) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 3, alpha = 0.5) +
  geom_text(
    data = top10_7day,
    aes(label = paste0("'", substr(year, 3, 4))),
    hjust = -0.5, vjust = 0.5, size = 2.5
  ) +
  geom_hline(
    data = rp_7day,
    aes(yintercept = roll7_sum),
    color = "#D62828", linetype = "dashed", linewidth = 1
  ) +
  geom_label(
    data = rp_7day,
    aes(x = Inf, y = roll7_sum,
        label = paste0(period_label, "\n~", round(rp_years, 1), "-yr")),
    hjust = 1.05, vjust = 0.5, size = 2.5,
    label.size = 0, fill = "white", color = "#D62828", fontface = "bold"
  ) +
  scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50)) +
  facet_wrap(~adm1_name, scales = "free_y") +
  labs(
    title = "Return Period: 7-Day Sum",
    subtitle = "January 2026 event in red",
    x = "Return Period (years)",
    y = "7-Day Sum (mm)"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold", size = 11),
    panel.grid.minor = element_blank()
  )

Summary

Takeaway: Table summarizes return period estimates by province. Higher return periods = rarer events.

summary_table <- bind_rows(
  rp_3day |>
    mutate(Metric = "3-Day") |>
    select(Province = adm1_name, Metric, Period = period_label,
           `Rainfall (mm)` = roll3_sum, `Return Period (yr)` = rp_years),
  rp_7day |>
    mutate(Metric = "7-Day") |>
    select(Province = adm1_name, Metric, Period = period_label,
           `Rainfall (mm)` = roll7_sum, `Return Period (yr)` = rp_years)
) |>
  mutate(`Rainfall (mm)` = round(`Rainfall (mm)`, 1),
         `Return Period (yr)` = round(`Return Period (yr)`, 1)) |>
  arrange(Province, Metric)

knitr::kable(summary_table)

Province	Metric	Period	Rainfall (mm)	Return Period (yr)
Gaza	3-Day	Jan 11-13	161.7	5.0
Gaza	7-Day	Jan 10-16	268.9	10.0
Maputo	3-Day	Jan 13-15	106.8	2.5
Maputo	7-Day	Jan 12-18	207.0	5.0
Sofala	3-Day	Jan 09-11	115.2	1.5
Sofala	7-Day	Jan 11-17	242.9	5.0

Flood Exposure Analysis (FloodScan)

df_lookup_adm2 <- df_lookup |> 
  filter(
    ADM1_PCODE %in% df_lookup_aoi$ADM1_PCODE,
    ADM_LEVEL ==2
  )

floodscan_adm2_unlabelled <- tbl(con, in_schema("app", "floodscan_exposure")) |>                           
    filter(pcode %in% df_lookup_adm2$ADM2_PCODE) |>                             
    collect()

floodscan_adm2 <- floodscan_adm2_unlabelled |> 
  left_join(
    df_lookup_adm2 |> 
      select(adm1_name=ADM1_NAME, adm2_name = ADM2_NAME, pcode = "ADM2_PCODE")
  )

FloodScan Exposure by District

Method: FloodScan flood extent crossed with WorldPop population to estimate people exposed per district.

Key finding: Heatmap shows temporal evolution of flood exposure. Districts ordered by total exposure - worst-affected at top. Darkest colors = highest population exposure.

Explore population exposure data interactively here

floodscan_recent <- floodscan_adm2 |>
  filter(valid_date >= max(valid_date) - 30) |>
  complete(
    valid_date = seq.Date(min(valid_date), max(valid_date), by = "day"),
    nesting(adm1_name, adm2_name),
    fill = list(sum = NA)
  )

ggplot(floodscan_recent, aes(x = valid_date, y = reorder(adm2_name, sum, \(x) sum(x, na.rm = TRUE)), fill = sum)) +
  geom_tile(color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(
    option = "inferno",
    direction = -1,
    labels = scales::comma,
    name = "Flood Exposure\n(sum)"
  ) +
  scale_x_date(date_breaks = "1 day", date_labels = "%b %d", expand = c(0, 0)) +
  facet_wrap(~adm1_name, ncol = 1, scales = "free_y") +
  labs(
    title = "FloodScan Exposure by District",
    subtitle = "Last 30 days | Districts ordered by total exposure",
    x = NULL,
    y = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(
    strip.text = element_text(face = "bold", size = 12),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 9),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 7),
    legend.position = "right",
    legend.key.height = unit(1.5, "cm")
  )

FloodScan Exposure Table

Takeaway: Same data in tabular form with exact values. Color intensity indicates exposure magnitude.

library(gt)

floodscan_wide <- floodscan_recent |>
  mutate(date_label = format(valid_date, "%b %d")) |>
  select(adm1_name, adm2_name, date_label, sum) |>
  pivot_wider(
    names_from = date_label,
    values_from = sum
  ) |>
  arrange(adm1_name, adm2_name)

date_cols <- floodscan_wide |>
  select(-adm1_name, -adm2_name) |>
  names()

sum_range <- range(floodscan_recent$sum, na.rm = TRUE)

floodscan_wide |>
  gt(groupname_col = "adm1_name") |>
  fmt_number(
    columns = all_of(date_cols),
    decimals = 0,
    use_seps = TRUE
  ) |>
  sub_missing(missing_text = "-") |>
  data_color(
    columns = all_of(date_cols),
    palette = "inferno",
    reverse = TRUE,
    na_color = "white",
    domain = sum_range
  ) |>
  cols_label(adm2_name = "District") |>
  tab_header(
    title = "FloodScan Exposure by District",
    subtitle = "Daily values (last 30 days)"
  ) |>
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_body()
  ) |>
  tab_style(
    style = cell_text(size = px(9)),
    locations = cells_column_labels()
  ) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_row_groups()
  ) |>
  tab_options(
    table.font.size = px(11),
    row_group.background.color = "#f0f0f0",
    column_labels.font.weight = "bold"
  )

FloodScan Exposure by District
Daily values (last 30 days)
District	Dec 18	Dec 19	Dec 20	Dec 21	Dec 22	Dec 23	Dec 24	Dec 25	Dec 26	Dec 27	Dec 28	Dec 29	Dec 30	Dec 31	Jan 01	Jan 02	Jan 03	Jan 04	Jan 05	Jan 06	Jan 07	Jan 08	Jan 09	Jan 10	Jan 11	Jan 12	Jan 13	Jan 14	Jan 15	Jan 16	Jan 17
Gaza
Bilene	18	15	16	17	16	16	18	19	17	17	36	15	-	-	19	19	18	40	40	43	27	26	48	-	14	4	20	43	579	917	2,767
Chibuto	436	890	400	402	438	505	517	495	479	468	469	364	-	-	1,546	1,545	1,528	1,635	1,619	1,707	1,552	1,524	1,525	-	560	558	1,720	1,561	1,738	2,638	33,785
Chicualacuala	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	148	140	57	0
Chigubo	0	0	0	0	0	0	0	2	10	8	7	1	-	-	0	0	0	0	0	0	0	0	0	-	21	0	6	103	74	85	346
Chokwe	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	9,640	16,666	14,825	18,760	46,454
Chongoene	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	3,667
Cidade De Xai-Xai	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	263
Guija	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	6	6	0	0	0	-	0	0	2,050	6,039	7,194	10,016	14,028
Limpopo	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	13,679
Mabalane	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	747	555	1,715	1,916	2,125	936
Mandlakaze	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	0
Mapai	0	0	0	0	0	0	0	12	44	34	33	10	-	-	0	0	0	0	0	0	0	0	0	-	0	0	2	459	303	339	346
Massangena	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	36	62	45	33	33
Massingir	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	1,451	1,890	2,303	1,150
Maputo
Boane	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	2,199	10,571
Cidade Da Matola	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	19,773	324
Magude	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	1,541	2,200	2,921	4,788
Manhiça	1,970	2,687	3,130	3,429	3,667	6,188	6,550	3,741	3,668	3,560	3,529	3,350	-	-	3,999	3,669	6,709	3,926	3,731	4,085	5,408	5,559	6,038	-	2,833	2,398	3,334	19,167	21,240	22,315	39,715
Marracuene	0	0	0	0	0	26	29	0	0	0	0	0	-	-	0	0	68	0	0	0	23	29	22	-	0	0	0	75	73	29	29
Matutuine	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	154	184	92	66	114	-	0	0	0	0	0	0	910
Moamba	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	938
Namaacha	0	0	0	0	0	0	0	0	0	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	551
Sofala
Buzi	0	0	4,491	5,308	6,216	12,184	12,592	13,600	13,608	27,536	29,191	32,778	-	-	25,641	23,303	16,597	14,969	13,463	12,741	12,178	11,758	11,920	-	12,815	15,456	19,234	24,266	29,342	27,804	29,733
Caia	0	0	0	0	4,312	4,690	4,619	1,304	1,363	10,779	12,626	12,538	-	-	10,842	9,911	7,547	6,767	6,751	6,728	2,136	1,914	4,907	-	5,850	5,761	5,496	5,756	5,549	4,823	5,099
Chemba	0	0	0	0	0	0	0	0	0	0	0	300	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	0
Cheringoma	0	0	0	0	80	94	102	123	102	386	418	323	-	-	188	212	166	183	197	242	205	181	182	-	240	232	136	140	197	175	256
Chibabava	0	0	0	0	0	0	188	345	1,010	1,153	1,105	225	-	-	0	0	0	0	0	0	247	207	0	-	692	4,002	3,213	3,066	2,915	2,206	2,398
Cidade Da Beira	0	0	11,839	40,889	72,583	43,471	34,329	24,072	1,567	37,426	39,272	33,004	-	-	46,434	33,898	12,659	2,203	10,906	2,164	10,199	1,938	5,897	-	19,880	15,773	17,474	2,575	7,204	1,244	17,393
Dondo	103	3,911	7,150	2,815	1,201	4,228	6,912	10,427	4,628	5,972	11,296	8,714	-	-	10,548	9,557	4,970	6,862	4,765	4,068	4,297	3,776	4,525	-	4,502	4,740	4,425	3,649	4,691	4,301	5,218
Gorongosa	0	237	791	837	573	835	1,387	2,734	2,483	4,759	4,755	3,798	-	-	2,934	2,586	2,269	1,985	2,043	2,108	2,138	2,144	1,715	-	2,136	2,428	4,634	2,892	2,945	2,529	2,926
Machanga	0	0	0	0	0	0	0	0	53	0	0	0	-	-	0	0	0	0	0	0	0	0	0	-	465	939	787	596	509	392	466
Maringue	0	0	0	0	0	0	0	839	975	31	30	0	-	-	0	0	0	0	0	0	0	0	0	-	0	0	0	0	0	0	0
Marromeu	0	0	0	0	93	109	118	137	112	2,555	5,861	7,521	-	-	208	213	483	207	218	264	202	175	179	-	157	143	148	154	160	140	142
Muanza	0	0	43	43	38	52	115	98	135	655	875	879	-	-	879	775	704	637	666	697	732	733	548	-	614	705	710	711	833	820	901
Nhamatanda	265	3,571	6,293	6,648	6,159	5,355	7,946	12,260	14,011	28,662	30,735	31,207	-	-	26,718	24,455	20,420	19,053	17,820	17,375	15,620	15,454	13,188	-	15,151	16,703	17,799	17,290	19,356	19,175	23,561

Flood Exposure Maps

Method: District-level population exposure mapped for the last 3 days. Gaza, Maputo, and Sofala provinces outlined in red.

Key finding: Top 3 most-affected districts labeled per province. Exposure concentrated in southern Mozambique, particularly along river valleys and low-lying coastal areas.

library(sf)
library(ggrepel)

sf_moz <- cumulus::download_fieldmaps_sf(iso3 = "moz", layer = "moz_adm2")

# Last 3 days of exposure data
floodscan_last3 <- floodscan_adm2 |>
  filter(valid_date >= max(valid_date) - 2)

# Join exposure to spatial (all provinces), filter out NA dates
sf_exposure <- sf_moz$moz_adm2 |>
  left_join(floodscan_last3, by = c("ADM2_PCODE" = "pcode")) |>
  filter(!is.na(valid_date))

# All province boundaries (light)
sf_provinces <- sf_moz$moz_adm2 |>
  group_by(ADM1_PT) |>
  summarise(geom = st_union(geom))

# Highlight boundaries for worst-affected provinces
sf_highlight <- sf_provinces |>
  filter(ADM1_PT %in% c("Gaza", "Maputo", "Sofala"))

# Province label positions (centroids)
sf_prov_labels <- sf_highlight |>
  mutate(
    centroid = st_centroid(geom),
    x = st_coordinates(centroid)[, 1],
    y = st_coordinates(centroid)[, 2],
    label = toupper(ADM1_PT)
  )

# Top 3 worst districts per province (for labeling)
sf_top3 <- sf_exposure |>
  filter(ADM1_PT %in% c("Gaza", "Maputo", "Sofala")) |>
  group_by(valid_date, ADM1_PT) |>
  slice_max(sum, n = 3) |>
  ungroup() |>

  mutate(centroid = st_centroid(geom)) |>
  mutate(
    x = st_coordinates(centroid)[, 1],
    y = st_coordinates(centroid)[, 2]
  )

ggplot() +
  geom_sf(data = sf_exposure, aes(fill = sum), color = "white", linewidth = 0.3) +
  geom_sf(data = sf_provinces, fill = NA, color = "grey70", linewidth = 0.5) +
  geom_sf(data = sf_highlight, fill = NA, color = "#D62828", linewidth = 1.5) +
  geom_text(
    data = sf_prov_labels,
    aes(x = x, y = y, label = label),
    size = 5, color = "grey50", fontface = "bold", alpha = 0.7
  ) +
  geom_label_repel(
    data = sf_top3,
    aes(x = x, y = y, label = ADM2_PT),
    size = 4,
    fill = alpha("white", 0.7),
    label.size = 0.3,
    fontface = "bold",
    segment.color = "grey40",
    segment.size = 0.4,
    box.padding = 0.5,
    point.padding = 0.3,
    max.overlaps = 20,
    min.segment.length = 0
  ) +
  scale_fill_viridis_c(
    option = "inferno",
    direction = -1,
    labels = scales::comma,
    na.value = "grey90",
    name = "Population Exposed"
  ) +
  facet_wrap(~valid_date, ncol = 1) +
  labs(
    title = "FloodScan Exposure by District",
    subtitle = "Gaza, Maputo, and Sofala highlighted | Top 3 districts labeled per province",
    caption = "Data: FloodScan & WorldPop"
  ) +
  theme_void(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold", size = 18, margin = margin(b = 10)),
    legend.position = "bottom",
    legend.key.width = unit(4, "cm"),
    legend.key.height = unit(0.6, "cm"),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 14, face = "bold"),
    plot.title = element_text(face = "bold", size = 22),
    plot.subtitle = element_text(color = "grey40", size = 14, margin = margin(b = 20))
  )