---
title: "How to use climate data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{How to use climate data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    # >>> Graphic settings
    # from https://r4ds.had.co.nz/graphics-for-communication.html?q=fig.asp#figure-sizing
    fig.width = 6,
    fig.asp = 0.618,
    fig.height= 6 * 0.618, # Needed for leaflet widget
    out.width = "100%",
    fig.align = "center"
)
```

```{r setup}
library(RADIS)
library(ggplot2)
library(tmap)

if (knitr::is_html_output()) {
  tmap_mode("view")
} else {
  tmap_mode("plot")
}
```

This vignette demonstrates how to extract and visualize climate data using the
RADIS package.
We extract daily precipitation and temperature from the SIM2 database on an area
covering municipalities around the city of Montpellier, France and compute mean
daily value for each municipality.

## Study area: municipalities around Montpellier, France

Let's define a bounding box roughly around Montpellier in RGF93 coordinate system:

```{r}
bbox_L93 <- sf::st_bbox(
  c(
    xmin = 760000,
    xmax = 783000,
    ymin = 6270000,
    ymax = 6290000
  ),
  crs = sf::st_crs(2154) # RGF93 / Lambert-93
)
```

And extract municipalities borders from the IGN API:

```{r}
rect <- sf::st_as_sfc(bbox_L93)
sf_M3M <- happign::get_wfs(
  x = rect,
  layer = "LIMITES_ADMINISTRATIVES_EXPRESS.LATEST:commune"
)
sf_M3M
```

The `sf` object `sf_M3M` is in WGS84 coordinate system, which is correct for
plotting it with leaflet:

```{r}
tm_basemap("OpenStreetMap.France") +
  tm_shape(sf_M3M) +
  tm_borders(col = "black", lwd = 1)
```

## Extracting raw climate data on the area

We download precipitation and mean temperature daily data from the SIM2 (SAFRAN)
database, which provides historical climate data for France using RADIS.api
(See https://api.g-eau.fr/__docs__).

The resulting variable is a `stars` object.

```{r}
s_climate <- get_sim2_daily(sf_M3M,
                               DATE__greater = as.Date("2025-06-15"),
                               DATE__less = as.Date("2025-07-15"),
                               fields = c("PRELIQ_Q", "T_Q"),
                               api_format = "ncdf",
                               overlay_mode = "extent")
s_climate
```

Let's see how the SAFRAN cells are intersecting the selected municipalities with
temperatures of the 15<sup>th</sup> day of the time series:

```{r}
#| map_climate_intersect,
#| fig.cap = "Climate data intersecting the municipalities around Montpellier"
plot_climate_intersect <- function(s, variable, i_time, limits = NULL) {
  if (get_spatial_format(s) == "raster") {
    data <- s[variable, , , i_time]
  } else {
    data <- s[variable, , i_time]
  }
  date <- stars::st_get_dimension_values(s, "time")[i_time]
  # Remove units which trigger ggplot2 Error in Ops.units(x, range[1])...
  for (i in seq_along(data)) {
    if (inherits(data[[i]], "units"))
      data[[i]] = units::drop_units(data[[i]])
  }

  ggplot() +
    geom_sf(data = sf::st_as_sf(data, long = TRUE),
            color = "grey",
            mapping = aes(fill = !!rlang::sym(names(data)[1]))) +
    geom_sf(
      data = sf_M3M,
      fill = NA,
      color = "grey20",
      linewidth = 1
    ) +
    coord_equal() +
    coord_sf(datum = sf::st_crs(s)) +
    geom_sf_text(data = sf_M3M, aes(label = nom_officiel), size = 2) +
    scale_fill_gradientn(colours = colorspace::diverge_hcl(7), limits = limits) +
    labs(title = paste0("Climate data: ", variable, " on ", date),
         fill = variable)
}
iTime <- 15
d <- s_climate$T_Q[, , iTime]
limits <- c(min(d, na.rm = TRUE), max(d, na.rm = TRUE))
plot_climate_intersect(s_climate, "T_Q", iTime, limits = limits)
```

## Climate data aggregation at municipalities scale

The aggregation can be operated with the function `spatial_combination` or
directly during data extraction with the argument `overlay_mode = "aggregate"` on
`get_sim2_daily` call.

```{r}
s_climate_agg <- spatial_combination(
  s_climate,
  sf_M3M,
  overlay_mode = "aggregate"
)
s_climate_agg
```

We can now show the aggregated temperatures at municipality scale:

```{r}
#| map_climate_agg,
#| fig.cap = "Climate data aggregated at municipalities scale"
plot_climate_intersect(s_climate_agg, "T_Q", iTime, limits = limits)
```

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