---
title: "v01-Get Started"
output:
  rmarkdown::html_vignette:
    css: styles/vignette.css
vignette: >
  %\VignetteIndexEntry{v01-Get Started}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup}
library(RADIS)
library(dplyr)
```

## Introduction

This vignette presents an example of how RADIS can be used.

## Inputs

Mandatory inputs for the package are:

- an `sf` object representing the study area.
- the year of interest, to retrieve crop data.
- the column name in the `sf` object that contains the unique identifier for
spatial features.
- the period of interest, to retrieve weather data.

## Study area

The study area used for demonstration and testing is provided as a shapefile
(`test.shp`) included in the package’s `extdata` directory.
It is loaded as an `sf` object and serves as the spatial reference for
extracting soil and climate data.

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

## Retrieve rpg data

The figure below show 2020 RPG data extracted for each polygon in the study
area.
This information is used to get crop types, plot areas, etc.

```{r rpg_data}
rpg <- get_rpg_data(sf = study_area_sf, year = 2020, id = "id")
plot(rpg[, "code_cultu"])
```

## Retrieve soil depth

The table below shows the soil depth values extracted for each spatial unit in
the study area.
These values are used as a basis for estimating available water capacity and
simulating water balance.

```{r soil_depth}
soil_depth <- get_soil_depth(sf = study_area_sf) %>%
  relocate(id)
knitr::kable(soil_depth, format = "pipe")
```

## Retrieve Soil Available Water Capacity (AWC)

### From BDGSF

AWC values from the BDGSF database are extracted for each polygon in the study
area.
This information is used to characterize the soil available water capacity
according to the national typology.

```{r soil_awc_bdgsf}
soil_awc_bdgsf <- get_soil_awc(sf = study_area_sf, source = "BDGSF") %>%
  relocate(id)
knitr::kable(soil_awc_bdgsf, format = "pipe")
```

### From Info&Sols

To use the Info&Sols database, you need to:

- Create an account in <https://entrepot.recherche.data.gouv.fr>.
- Get an API token from this website (click on your name in the upper right
corner and go to `API Token`).
- Set an environment variable named `API_TOKEN_DATAVERSE` or pass you API token
directly as the `key` argument in the `get_soil_awc` function.

RADIS can retrieve soil AWC :

- With coarse elements:

AWC estimates from the Info&Sols database, including the effect of coarse
elements, provide a more realistic value of available water capacity,
especially in stony soils.

```{r soil_awc_infosols_with_coarse}
soil_awc_infosols_with_coarse <- get_soil_awc(
  sf = study_area_sf, source = "Info&Sols", key = ""
)
knitr::kable(soil_awc_infosols_with_coarse, format = "pipe")
```

- Without coarse elements

AWC values extracted without accounting for coarse elements represent a
theoretical maximum, useful for comparing the impact of coarse fragments on
water retention.

```{r soil_awc_infosols_without_coarse}
soil_awc_infosols_without_coarse <- get_soil_awc(
  sf = study_area_sf, source = "Info&Sols", with_coarse_elements = FALSE
)
knitr::kable(soil_awc_infosols_without_coarse, format = "pipe")
```

## Retrieve Soil Texture

The `get_soil_texture()` function extracts soil texture information for each
spatial unit in the study area.
This function supports multiple data sources (e.g., Info&Sols, HWSD, SoilGrids)
 and allows you to specify the overlay mode for aggregating values.
The resulting data includes proportions of sand, silt, and clay, which are essential for characterizing soil physical properties.

### From Info&Sols

```{r soil_texture_infosols}
# Raw format
soil_texture_infosols <- get_soil_texture(
  study_area_sf,
  source = "infosols",
  overlay_mode = "aggregate"
)
knitr::kable(soil_texture_infosols, format = "pipe")
```

### From HWSD

```{r soil_texture_hwsd}
# Raw format
soil_texture_hwsd <- get_soil_texture(
  study_area_sf,
  source = "hwsd",
  overlay_mode = "aggregate"
)
knitr::kable(soil_texture_hwsd, format = "pipe")
```

# SoilGrids

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

This is a comparison of soil texture results obtained from different data
sources (Info&Sols, HWSD, SoilGrids).
This helps to highlight differences in estimates and potential implications for
downstream analyses.

- **Info&Sols**: National French database, high spatial resolution, based on
field observations and expert knowledge.
- **HWSD**: Global database, moderate resolution, harmonized from multiple
sources.
- **SoilGrids**: Global database, high spatial resolution, based on machine
learning predictions.

Comparing these sources allows to assess the variability and reliability of
soil texture data for their specific study area.

```{r plot_soil_texture}
# 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")
  )
```

## Extracted Data Vizualisation

```{r map, eval=TRUE, fig.width=15}
merged <- Reduce(
  f = function(b, lyr) {
    attrs <- sf::st_set_geometry(lyr, NULL)
    left_join(b, attrs, by = "id")
  },
  x = list(
    soil_depth,
    soil_awc_bdgsf,
    soil_awc_infosols_with_coarse,
    soil_awc_infosols_without_coarse,
    soil_texture_infosols
  ),
  init = study_area_sf
)

view_data_map(
  sf = merged,
  cols_to_show = c(
    "AREA.x", "soil_depth", "argile.0_5", "limon.0_5", "sable.0_5",
    "awc_mm_015_030", "awc_mm_030_060", "awc_mm_060_100"
  ),
  rpg_sf = rpg,
  basemaps = "OpenStreetMap"
)
```

```{r, include=FALSE}
# Clean up cache
clean_radis_cache()
```