study_area_sf <- sf::read_sf(
  system.file("extdata/study_area/test.shp", package = "RADIS")
)
plot(study_area_sf)

rpg <- get_rpg_data(sf = study_area_sf, year = 2020, id = "id")

plot(rpg[, "code_cultu"])

soil_depth <- get_soil_depth(sf = study_area_sf) %>%
  relocate(id)

knitr::kable(soil_depth, format = "pipe")

# Raw format
soil_texture_soilgrids <- get_soil_texture(
  study_area_sf,
  source = "soilgrids",
  overlay_mode = "aggregate"
)
knitr::kable(soil_texture_soilgrids, format = "pipe")

# reshape InfoSols
tb_infosols <- tibble::tibble(
  depth = c("0_5", "5_15", "15_30", "30_60", "60_100", "100_200"),
  top   = c(0, 5, 15, 30, 60, 100),
  bot   = c(5, 15, 30, 60, 100, 200)
)
df_infosols <- soil_texture_infosols %>%
  tidyr::pivot_longer(
    matches("^(argile|sable|limon)\\.\\d+_\\d+$"),
    names_to  = c("texture", "depth"),
    names_sep = "\\.",
    values_to = "pc"
  ) %>%
  dplyr::mutate(texture = dplyr::recode(
    texture,
    argile = "Clay",
    limon = "Silt",
    sable = "Sand"
  )) %>%
  dplyr::left_join(tb_infosols, by = "depth") %>%
  dplyr::mutate(source = "InfoSols")

# reshape HWSD
tb_hwsd <- tibble::tibble(
  d   = paste0("D", 1:7),
  top = c(0, 20, 40, 60, 80, 100, 150),
  bot = c(20, 40, 60, 80, 100, 150, 200)
)
df_hwsd <- soil_texture_hwsd %>%
  tidyr::pivot_longer(
    matches("^(SAND|SILT|CLAY)_D\\d+$"),
    names_to      = c("texture", "d"),
    names_pattern = "^(.*)_(D\\d+)$",
    values_to     = "pc"
  ) %>%
  dplyr::mutate(texture = dplyr::recode(
    texture,
    CLAY = "Clay",
    SILT = "Silt",
    SAND = "Sand"
  )) %>%
  dplyr::left_join(tb_hwsd, by = "d") %>%
  dplyr::mutate(source = "HWSD")

# reshape SoilGrids
tb_soilgrids <- tibble::tibble(
  depth = c("0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm"),
  top   = c(0, 5, 15, 30, 60, 100),
  bot   = c(5, 15, 30, 60, 100, 200)
)
df_soilgrids <- soil_texture_soilgrids %>%
  tidyr::pivot_longer(
    matches("^(sand|silt|clay)_\\d+-\\d+cm_mean$"),
    names_to      = c("texture", "depth"),
    names_pattern = "^(sand|silt|clay)_(\\d+-\\d+cm)_mean$",
    values_to     = "pc"
  ) %>%
  dplyr::mutate(texture = dplyr::recode(
    toupper(texture),
    CLAY = "Clay",
    SILT = "Silt",
    SAND = "Sand"
  )) %>%
  dplyr::left_join(tb_soilgrids, by = "depth") %>%
  dplyr::mutate(source = "SoilGrids")


# combine + prep
df <- dplyr::bind_rows(df_infosols, df_hwsd, df_soilgrids) %>%
  dplyr::filter(!is.na(pc), !is.na(top), !is.na(bot)) %>%
  dplyr::group_by(source, texture, top, bot) %>%
  dplyr::summarise(pc = mean(pc, na.rm = TRUE), .groups = "drop") %>%
  dplyr::mutate(
    y_mid  = (top + bot) / 2,
    source = factor(source, c("InfoSols", "HWSD", "SoilGrids")),
    x_num  = as.numeric(source)
  )

# plot
ggplot2::ggplot(df, ggplot2::aes(x = source)) +
  ggplot2::geom_rect(
    ggplot2::aes(
      xmin = x_num - 0.35, xmax = x_num + 0.35,
      ymin = top, ymax = bot, fill = pc
    ),
    color = "grey30", linewidth = 0.2
  ) +
  ggplot2::geom_text(
    ggplot2::aes(y = y_mid, label = paste0(round(pc, 1), "%")),
    color = "black", size = 3
  ) +
  ggplot2::facet_wrap(~texture, nrow = 1) +
  ggplot2::scale_y_reverse(
    name = "Depth (cm)",
    limits = c(max(df$bot), 0),
    expand = c(0, 0)
  ) +
  ggplot2::scale_x_discrete(name = "Source") +
  ggplot2::scale_fill_gradient(
    low = "#deebf7",
    high = "#3182bd",
    name = "Percentage"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    panel.grid.major.x = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    strip.text = ggplot2::element_text(face = "bold")
  )