R st_nearest_point returning linestring not touching second geometry

I'm using the st_nearest_point function to find the closest spot on a line to a nearby point. I am not sure what the issue is, but I cannot get the function to return a linestring with any point touching the line. I need to use the 4326 crs because that is what leaflet maps use and in my real application my starting point is a map_click capture. st_nearest_point returns a linestring that is close but not touching the line. What gives? Please see example below.

library(sf)
library(dplyr)
library(mapview)

#Creating some dummy data to test with
# Create a line object (LINESTRING)
line <- st_sfc(st_linestring(matrix(c(-122.4, 37.7, -122.5, 37.8), ncol = 2, byrow = TRUE)), crs = 4326)
# Create a point object (POINT) that is close to but does not touch the line
point <- st_sfc(st_point(c(-122.45, 37.78)), crs = 4326)
# Create sf objects
line_sf <- st_sf(geometry = line)
point_sf <- st_sf(geometry = point)


#Plotting the objects so you can see they don't touch
plot(st_geometry(line_sf), col = 'blue', lwd = 2)
plot(st_geometry(point_sf), col = 'red', add = TRUE)

#Confirming with st_touches
st_touches(line_sf,point_sf)


#Lets find the closest point on the line and use that

geo <- st_nearest_points(point_sf, line_sf)
lngs = geo |> st_length() # Next few lines are redundant here but important if there are multiple lines
nearestp = min(lngs)
my_linestring <- geo[nearestp == lngs]
new_point <- st_cast(my_linestring, 'POINT')[c(FALSE,TRUE)] |> st_as_sf()


#Now lets plot this new point   
plot(st_geometry(line_sf), col = 'blue', lwd = 2)
plot(st_geometry(point_sf), col = 'red', add = TRUE)
plot(st_geometry(new_point), col = 'green', add = TRUE)

#Success? Not so fast lets confirm with st_touches

st_touches(line_sf,new_point)

#Lets use mapview to see if we can see by zooming in the farthest we can

mapview(line_sf) + mapview(new_point, col.region  = "green") + mapview(point_sf, col.region  = "red")

#Looks good but if you zoom in you will see the points indeed do not touch. 

What gives?

This is probably a shortcoming of floating point arithmetic as discussed here: https://github.com/r-spatial/sf/issues/1503

There is a workaround using sfnetworks https://luukvdmeer.github.io/sfnetworks/articles/sfn03_join_filter.html#blending-points-into-a-network

That's the approach taken in the code below. While the point does end up on the line, st_touches() still returns false or empty.

library(sfnetworks)

#blend the line & points into an sf_newtork
blended <- st_network_blend(as_sfnetwork(line_sf), new_point)
#select the edges & return a linestring
line_n <- blended %>% activate('edges') %>% st_as_sf() %>% st_union() %>% st_as_sf() %>% st_cast('LINESTRING')
# select the nodes & return as points
points_n <- blended %>% activate("nodes") %>% st_as_sf() %>% st_combine() %>% st_as_sf() %>% st_cast('POINT')
point_n <- points_n[3,] #just picking the right point (not the end nodes)

mapview(line_n) + mapview(points_n, col.region = 'green') + mapview(point_n, col.region = 'red')

st_touches(line_n, point_n)
#Sparse geometry binary predicate list of length 1, where the predicate was #`touches'
# 1: (empty)