Skip to content

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")