library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(lwgeom)
## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 7.0.0
library(nngeo)
library(ggplot2)

# Sample data
url = "https://geobgu.xyz/web-mapping2/additional/flow_direction/lines.geojson"
dat = st_read(url)
## Reading layer `lines2' from data source `https://geobgu.xyz/web-mapping2/additional/flow_direction/lines.geojson' using driver `GeoJSON'
## Simple feature collection with 102 features and 6 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 0.00035 ymin: 0.00025 xmax: 0.03225 ymax: 0.02105
## geographic CRS: WGS 84
dat = dat[, "geometry"]
st_crs(dat) = NA
elev = st_make_grid(dat, n = c(30, 30))
elev = st_as_sf(elev)
elev$elev = st_coordinates(st_centroid(elev))[,1]

# Plot
plot(elev, reset = FALSE)
plot(st_geometry(dat), add = TRUE)

# Extract elevation
dat$from_elev = NA
dat$to_elev = NA
for(i in 1:nrow(dat)) {
    start = st_startpoint(dat[i, ])
    dat$from_elev[i] = elev[start, ]$elev
    end = st_endpoint(dat[i, ])
    dat$to_elev[i] = elev[end, ]$elev
}

# Set direction
sel = dat$from_elev > dat$to_elev
dat = rbind(
    dat[!sel, ],
    st_reverse(dat[sel, ])
)

# To segments
seg = st_segments(dat, progress = FALSE)

# To coordinates
start = st_coordinates(st_startpoint(seg))
end = st_coordinates(st_endpoint(seg))
coords = cbind(start, end)
coords = as.data.frame(coords)
names(coords) = c("from_x", "from_y", "to_x", "to_y")
start = st_coordinates(st_startpoint(dat))
end = st_coordinates(st_endpoint(dat))
coords2 = cbind(start, end)
coords2 = as.data.frame(coords2)
names(coords2) = c("from_x", "from_y", "to_x", "to_y")

# Plot
ggplot() +
    geom_segment(data = coords, aes(x = from_x, y = from_y, xend = to_x, yend = to_y), colour = "grey") +
    geom_segment(data = coords2, aes(x = from_x, y = from_y, xend = to_x, yend = to_y),arrow = arrow(length = unit(0.03, "npc")))