Improve bug report or checks for extract_reaches
When you use extract_reaches but have no subcatchment, then (logically) extract_reaches does not work, as there are no reaches.
We have the following error:
> reaches <- extract_reaches(sb = sb, overwrite = overwrite)
Reading layer `streams_55944a0b90f73bd089bcdd805f1d001a' from data source
`C:\Users\gthirel\AppData\Local\Temp\airGRccia_cache\streams_55944a0b90f73bd089bcdd805f1d001a.shp' using driver `ESRI Shapefile'
Simple feature collection with 20127 features and 1 field
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 714437 ymin: 6236544 xmax: 778416.6 ymax: 6356142
Projected CRS: RGF93 v1 / Lambert-93
Erreur dans sb$mvd_outlet_x[i_downs] : type 'list' d'indice incorrect
In my opinion, there should be a check in extract_reaches to verify that we have subcatchments. Maybe checking the dimension of sb and outputing a bug report stating that sb must contain more than one catchment would be the easiest and most comprehensible.
Reproducible code:
overwrite <- TRUE
bbox <- sf::st_bbox(
c(
xmin = 716000,
ymin = 6238000,
xmax = 774000,
ymax = 6355000
),
crs = 2154
)
# Create an sf object from the bounding box
bbox_sf <- sf::st_as_sfc(bbox) %>% st_as_sf()
tmap_mode("view")
#> ℹ tmap mode set to "view".
tm_basemap("OpenTopoMap") +
tm_shape(bbox_sf) +
tm_polygons(lwd = 3, col = "red", fill = NA, fill_alpha = 0)
dem_file <- getCachePath(bbox, "dem", "tif")
dem <- get_ign_dem(bbox, filename = dem_file, overwrite = overwrite)
file_streams <- getCachePath(bbox, "stream", "shp")
streams <- get_ign_streams(bbox, filename = file_streams, overwrite = overwrite)
original_mode <- tmap_mode()
tmap_mode("plot")
tm_basemap("OpenStreetMap.France") +
tm_shape(streams) +
tm_lines(col = "nature")
tmap_mode(original_mode)
stations <- hubeau::get_hydrometrie_sites(
code_site = c(
"Y2140010"
)
)
stations
file_sb <- getCachePath("V01", "sb", "RDS")
outlets <- stations[, c("code_site", "coordonnee_x_site", "coordonnee_y_site")]
colnames(outlets) <- c("id", "x", "y")
sb <- extract_sub_basins(outlets = outlets, dem = dem, streams = streams, filename = file_sb, overwrite = overwrite)
plot(terra::vect(sb), c("id", "down"))
reaches <- extract_reaches(sb = sb, overwrite = overwrite)
moved_outlets <- sf::st_read(attr(sb, "files")["moved_outlets"])
tm_basemap("OpenStreetMap.France") +
tm_shape(sb) +
tm_polygons(fill = "id", fill_alpha = 0.3) +
tm_shape(reaches) +
tm_lines(col = "blue", lwd = 3) +
tm_shape(moved_outlets) +
tm_symbols(fill = "id")