If flowlines aren't digitized in the expected direction, this will reorder the nodes so they are.
fix_flowdir(comid, network = NULL, fn_list = NULL)Arguments
- comid
integer The COMID of the flowline to check
- network
sf data.frame The entire network to check from. Requires a "toCOMID" field.
- fn_list
list containing named elements `flowline`, `network`, and `check_end`, where `flowline` is the flowline to be checked and `network` the feature up or downstream of the flowline to be checked, and `check_end` is `"start"` or `"end"` depending if the `network` input is upstream (`"start"`) or downstream (`"end"`) of the flowline to be checked. This option allows pre-compilation of pairs of features which may be useful for very large numbers of flow direction checks.
Value
a geometry for the feature that has been reversed if needed.
Examples
source(system.file("extdata/sample_data.R", package = "nhdplusTools"))
fline <- sf::read_sf(sample_data, "NHDFlowline_Network")
# We add a tocomid with prepare_nhdplus
fline <- sf::st_sf(prepare_nhdplus(fline, 0, 0, 0, FALSE),
geom = sf::st_zm(sf::st_geometry(fline)))
#> Warning: removing geometry
#> Warning: Got NHDPlus data without a Terminal catchment. Attempting to find it.
#> Warning: Removed 0 flowlines that don't apply.
#> Includes: Coastlines, non-dendritic paths,
#> and networks with drainage area less than 0 sqkm, and drainage basins smaller than 0
# Look at the end node of the 10th line.
(n1 <- get_node(fline[10, ], position = "end"))
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -89.34813 ymin: 42.99555 xmax: -89.34813 ymax: 42.99555
#> Geodetic CRS: GRS 1980(IUGG, 1980)
#> # A tibble: 1 × 1
#> geometry
#> * <POINT [°]>
#> 1 (-89.34813 42.99555)
# Break the geometry by reversing it.
sf::st_geometry(fline)[10] <- sf::st_reverse(sf::st_geometry(fline)[10])
# Note that the end node is different now.
(n2 <- get_node(fline[10, ], position = "end"))
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -89.36355 ymin: 43 xmax: -89.36355 ymax: 43
#> Geodetic CRS: GRS 1980(IUGG, 1980)
#> # A tibble: 1 × 1
#> geometry
#> * <POINT [°]>
#> 1 (-89.36355 43)
# Pass the broken geometry to fix_flowdir with the network for toCOMID
sf::st_geometry(fline)[10] <- fix_flowdir(fline$COMID[10], fline)
# Note that the geometry is now in the right order.
(n3 <- get_node(fline[10, ], position = "end"))
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -89.34813 ymin: 42.99555 xmax: -89.34813 ymax: 42.99555
#> Geodetic CRS: GRS 1980(IUGG, 1980)
#> # A tibble: 1 × 1
#> geometry
#> * <POINT [°]>
#> 1 (-89.34813 42.99555)
plot(sf::st_geometry(fline)[10])
plot(n1, add = TRUE)
plot(n2, add = TRUE, col = "blue")
plot(n3, add = TRUE, cex = 2, col = "red")
