How to use climate data

library(RADIS)
library(ggplot2)
library(tmap)

if (knitr::is_html_output()) {
  tmap_mode("view")
} else {
  tmap_mode("plot")
}
#> ℹ tmap modes "plot" - "view"
#> ℹ toggle with `tmap::ttm()`

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:

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:

rect <- sf::st_as_sfc(bbox_L93)
sf_M3M <- happign::get_wfs(
  x = rect,
  layer = "LIMITES_ADMINISTRATIVES_EXPRESS.LATEST:commune"
)
sf_M3M
#> Simple feature collection with 43 features and 17 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 3.621409 ymin: 43.46662 xmax: 4.093463 ymax: 43.78872
#> Geodetic CRS:  WGS 84
#> # A tibble: 43 × 18
#>    id    cleabs nom_officiel nom_officiel_en_maju…¹ statut code_insee population
#>    <chr> <chr>  <chr>        <chr>                  <chr>  <chr>           <int>
#>  1 comm… COMMU… Argelliers   ARGELLIERS             Commu… 34012            1002
#>  2 comm… COMMU… Assas        ASSAS                  Commu… 34014            1405
#>  3 comm… COMMU… Baillargues  BAILLARGUES            Commu… 34022            7829
#>  4 comm… COMMU… Castelnau-l… CASTELNAU-LE-LEZ       Commu… 34057           26058
#>  5 comm… COMMU… Castries     CASTRIES               Commu… 34058            6883
#>  6 comm… COMMU… Clapiers     CLAPIERS               Commu… 34077            6006
#>  7 comm… COMMU… Le Crès      LE CRES                Commu… 34090            9226
#>  8 comm… COMMU… Combaillaux  COMBAILLAUX            Commu… 34082            2001
#>  9 comm… COMMU… Cournonterr… COURNONTERRAL          Commu… 34088            7359
#> 10 comm… COMMU… Fabrègues    FABREGUES              Commu… 34095            7447
#> # ℹ 33 more rows
#> # ℹ abbreviated name: ¹nom_officiel_en_majuscules
#> # ℹ 11 more variables: date_du_recensement <date>, organisme_recenseur <chr>,
#> #   code_insee_du_canton <chr>, code_insee_de_l_arrondissement <chr>,
#> #   code_insee_du_departement <chr>, code_insee_de_la_region <chr>,
#> #   codes_siren_des_epci <chr>, code_siren <chr>, code_postal <chr>,
#> #   superficie_cadastrale <int>, geometry <MULTIPOLYGON [°]>

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

tm_basemap("OpenStreetMap.France") +
  tm_shape(sf_M3M) +
  tm_borders(col = "black", lwd = 1)
#> Linking to GEOS 3.14.1, GDAL 3.12.2, PROJ 9.7.1; sf_use_s2() is FALSE

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.

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")
#> for reading time stamps, package CFtime is required: please install it first
#> Raster time series with transformed CRS results in a grid with curvilinear coordinates.
#> There is no satisfactory way to save this type of object in NetCDF format, use RDS instead.
s_climate
#> stars object with 3 dimensions and 2 attributes
#> attribute(s):
#>                 Min. 1st Qu. Median       Mean 3rd Qu. Max. NAs
#> PRELIQ_Q [mm/d]  0.0     0.0    0.0  0.7822581     0.1 38.7  62
#> T_Q [°C]        18.8    23.9   26.2 26.0905530    28.0 31.4  62
#> dimension(s):
#>      from to     offset  delta refsys point                values x/y
#> x       1  6         NA     NA WGS 84    NA [6x5] 3.571,...,4.073 [x]
#> y       1  5         NA     NA WGS 84    NA [6x5] 43.49,...,43.78 [y]
#> time    1 31 2025-06-15 1 days   Date FALSE                  NULL    
#> curvilinear grid

Let’s see how the SAFRAN cells are intersecting the selected municipalities with temperatures of the 15^th day of the time series:

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)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

Climate data intersecting the municipalities around Montpellier

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.

s_climate_agg <- spatial_combination(
  s_climate,
  sf_M3M,
  overlay_mode = "aggregate"
)
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar
s_climate_agg
#> stars object with 2 dimensions and 2 attributes
#> attribute(s):
#>                     Min. 1st Qu.   Median       Mean     3rd Qu.     Max.
#> PRELIQ_Q [mm/d]  0.00000     0.0  0.00000  0.6140612  0.09697435 34.45291
#> T_Q [°C]        20.15246    24.1 26.21598 26.1536181 28.05162634 30.49010
#> dimension(s):
#>          from to     offset  delta refsys point
#> geometry    1 43         NA     NA WGS 84 FALSE
#> time        1 31 2025-06-15 1 days   Date FALSE
#>                                                                 values
#> geometry MULTIPOLYGON (((3.622231 ...,...,MULTIPOLYGON (((3.825793 ...
#> time                                                              NULL

We can now show the aggregated temperatures at municipality scale:

plot_climate_intersect(s_climate_agg, "T_Q", iTime, limits = limits)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

Climate data aggregated at municipalities scale