V02 - Hydroclimatic data preprocessing
The final goal of this vignette is to create a BasinsObs object containing both hydrometric and climatic observations.
Retrieving and format climatic data
Climatic data: French SAFRAN reanalysis
First, we retrieve climatic data from the SAFRAN reanalysis produced by Meteo France (todo: add some details about this database), on the study area. We proceed by querying the RADIS API with the function RADIS package developed by the G-EAU Mixte Research Unit.
# Load sub-basins shapefile for getting the extent of the data retrieval
sf_sb <- readRDS(getCachePath("V01", "sb", "RDS"))
# Get the SAFRAN stars object with RADIS
s_safran <- get_climate_safran(sf_sb,
date_start = "2012-08-01",
date_end = "2024-12-31",
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.The s_safran object is a stars object
containing the climatic variables in grid format:
plot_climate_bvi <- function(s, sf_sb, variable = "TempMean", i_time, limits = c(-20, 40)) {
ggplot() +
geom_stars(data = s[variable, , , i_time]) +
geom_sf(
data = sf_sb,
fill = NA,
color = "grey20",
linewidth = 1
) +
#coord_equal() +
coord_sf(datum = sf::st_crs(s)) +
geom_sf_text(data = sf_sb, aes(label = id), size = 1.5) +
scale_fill_gradient2(
limits = limits,
low = scales::muted("blue"),
high = scales::muted("red"),
mid = "white"
) +
ggtitle(sprintf(
"Variable '%s' at date %s",
variable,
format(stars::st_get_dimension_values(s_safran, "time")[i_time])
))
}
# Every layer should be in same CRS in order to be plotted together
plot_climate_bvi(s_safran, sf_sb, i_time = 1)
Semi-distrbuted hydrologixcal models such as airGRiwrm require mean
climate data at sub-basin scale. Such data are available using the
agrument overlay_mode = "aggregate" when retrieving the
data:
s_sb <- get_climate_safran(
sf_sb,
date_start = "2012-08-01",
date_end = "2024-12-31",
overlay_mode = "aggregate"
)
#> for reading time stamps, package CFtime is required: please install it first
#> Warning in CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,
#> : GDAL Error 6: SetIndexingVariable() not implemented
#> Warning in CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,
#> : GDAL Error 6: SetIndexingVariable() not implemented
plot(s_sb["TempMean", , 1:6])
Creation of the dataset for hydrometric data
We retrieve observed flow time series from site codes of 8 character-lengths used in the sub-basin delineation.
Qobs <- get_hubeau_flows(
x = sf_sb$id, # Use sub-basin ids
date_start = "2012-08-01",
date_end = "2024-12-31"
)
#> Listing station codes ■■■■■■■■■■■ 33% | ETA: 2s
#> Listing station codes ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 100% | ETA: 0s
#>
#> Downloading station Y223001001 from 2012-08-01 to 2024-12-31
#> Warning in FUN(X[[i]], ...): No observation found for code_entite Y223001001
#> Downloading station Y221001001 from 2012-08-01 to 2024-12-31
#> Downloading station Y221001002 from 2012-08-01 to 2024-12-31
#> Warning in FUN(X[[i]], ...): No observation found for code_entite Y221001002
#> Downloading station Y210002001 from 2012-08-01 to 2024-12-31
#> Downloading station Y214002001 from 2012-08-01 to 2024-12-31
#> Downloading station Y214001001 from 2012-08-01 to 2024-12-31
#> Downloading station Y214001002 from 2012-08-01 to 2024-12-31
#> Downloading station Y230002001 from 2012-08-01 to 2024-12-31
summary(Qobs)
#> DatesR Y221001001 Y210002001
#> Min. :2012-08-01 00:00:00 Min. : 0.400 Min. : 1.111
#> 1st Qu.:2015-09-08 18:00:00 1st Qu.: 1.018 1st Qu.: 3.446
#> Median :2018-10-16 12:00:00 Median : 1.759 Median : 8.347
#> Mean :2018-10-16 12:00:00 Mean : 3.510 Mean : 17.448
#> 3rd Qu.:2021-11-23 06:00:00 3rd Qu.: 3.581 3rd Qu.: 17.241
#> Max. :2024-12-31 00:00:00 Max. :121.779 Max. :474.403
#> NAs :129
#> Y214002001 Y214001002 Y230002001 Y214001001
#> Min. : 0.692 Min. : 1.453 Min. : 1.874 Min. : 7.725
#> 1st Qu.: 3.792 1st Qu.: 3.774 1st Qu.: 6.047 1st Qu.: 13.579
#> Median : 10.312 Median : 11.333 Median : 13.964 Median : 22.386
#> Mean : 22.084 Mean : 24.211 Mean : 30.367 Mean : 42.079
#> 3rd Qu.: 23.477 3rd Qu.: 24.038 3rd Qu.: 27.853 3rd Qu.: 41.419
#> Max. :613.311 Max. :825.129 Max. :1011.325 Max. :804.745
#> NAs :475 NAs :11 NAs :3270We get the observed flows by station code (10 character-lengths) and there are potentially more than one station by site.
Nota Bene: no observed flows are available for the gauging site Y2230010 (Le Salagou à Clermont-l’Hérault [Mas Audran]).
We can plot the data availability for each station:
# Plot data availability for each station
ggplot(
Qobs %>% tidyr::pivot_longer(-DatesR, names_to = "code_station", values_to = "Q") %>% filter(!is.na(Q)),
aes(x = as.POSIXct(DatesR), y = code_station)
) +
geom_point() +
ggtitle("Gauging station data availabilities")
As station Y214001001 has definitely been replaced by
station Y214001002 in 2020 we remove the latter for this
site.
Then We can convert station code into site codes to be compliant with the code used with sub-basins.
Create the BasinsObs object
BasinsObs <- CreateBasinsObs(s_sb, ids = sf_sb$id, Qobs = Qobs)
str(BasinsObs)
#> List of 5
#> $ DatesR : POSIXlt[1:4536], format: "2012-08-01" "2012-08-02" ...
#> $ Precip :'data.frame': 4536 obs. of 6 variables:
#> ..$ Y2230010: num [1:4536] 0 0 0 0.324 3.24 ...
#> ..$ Y2210010: num [1:4536] 0 0 0 1.79 8.5 ...
#> ..$ Y2100020: num [1:4536] 0 0 0 6.67 13.53 ...
#> ..$ Y2140020: num [1:4536] 0 0 0 0.945 3.338 ...
#> ..$ Y2140010: num [1:4536] 0 0 0 0.474 3.602 ...
#> ..$ Y2300020: num [1:4536] 0 0 0 0.467 3.78 ...
#> $ PotEvap :'data.frame': 4536 obs. of 6 variables:
#> ..$ Y2230010: num [1:4536] 4.86 9 6.84 6.77 5.41 ...
#> ..$ Y2210010: num [1:4536] 5.76 8.43 7.04 6.76 5.71 ...
#> ..$ Y2100020: num [1:4536] 6.44 7.43 7.21 5.83 4.34 ...
#> ..$ Y2140020: num [1:4536] 6.24 7.1 6.27 6.21 5.29 ...
#> ..$ Y2140010: num [1:4536] 5.22 8.57 6.9 6.56 5.42 ...
#> ..$ Y2300020: num [1:4536] 4.99 8.84 7.18 6.66 5.43 ...
#> $ TempMean:'data.frame': 4536 obs. of 6 variables:
#> ..$ Y2230010: num [1:4536] 23.3 26 25.1 25.1 23 ...
#> ..$ Y2210010: num [1:4536] 21.5 23.1 22.2 22.2 20.3 ...
#> ..$ Y2100020: num [1:4536] 22.4 23.4 22.6 22.4 19.9 ...
#> ..$ Y2140020: num [1:4536] 23.1 25 25.5 24.8 23.1 ...
#> ..$ Y2140010: num [1:4536] 23.3 25.9 25.5 25.1 23.2 ...
#> ..$ Y2300020: num [1:4536] 23.4 26.3 25.6 25.3 23.3 ...
#> $ Qobs : tibble [4,536 × 5] (S3: tbl_df/tbl/data.frame)
#> ..$ Y2210010: num [1:4536] 0.769 0.666 0.661 0.71 0.702 0.726 0.722 0.697 0.7 0.669 ...
#> ..$ Y2100020: num [1:4536] 2.91 2.9 2.88 2.86 2.92 ...
#> ..$ Y2140020: num [1:4536] 1.61 1.55 1.55 1.5 1.52 ...
#> ..$ Y2140010: num [1:4536] 2.24 2.42 2.51 2.46 2.44 ...
#> ..$ Y2300020: num [1:4536] 3.67 3.67 3.56 3.39 3.4 ...
#> ..- attr(*, "codes_entite_station")='data.frame': 8 obs. of 2 variables:
#> .. ..$ code_entite : chr [1:8] "Y2230010" "Y2210010" "Y2210010" "Y2100020" ...
#> .. ..$ code_station: chr [1:8] "Y223001001" "Y221001001" "Y221001002" "Y210002001" ...
#> ..- attr(*, "DatesR")= POSIXlt[1:4536], format: "2012-08-01" "2012-08-02" ...
#> - attr(*, "class")= chr [1:2] "BasinsObs" "list"Save the BasinsObs object on the cloud