Basically the problem was found to be only at the value of 180.0, so by clamping to 179.9, all was solved and the maps display perfectly now. A short python script resolved this situation. It's amazing that after asking on the gdal-dev mailing list, other websites, this one, Twitter, and Facebook, not a single GIS expert was able to dispense this small piece of expertise smile.
If you snap the point to the line then you can split the line and then extract the parts from the resulting collection. Use a really small tolerance, I don't know how small it needs to be...
> site_snap = st_snap(site, reach, tol=1e-9)
> parts = st_collection_extract(st_split(reach$geometry, site_snap$geometry),"LINESTRING")
parts
is now two features, either side of your point:
> parts
Geometry set for 2 features
geometry type: LINESTRING
dimension: XY
bbox: xmin: 688358.8 ymin: 171876.4 xmax: 688465 ymax: 172171.2
epsg (SRID): NA
proj4string: +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs
LINESTRING (688400.1 171876.4, 688394.4 171890....
LINESTRING (688448.1 172086.7, 688445.9 172109....
Your original site
doesn't intersect the line, so your split returned nothing. Test with:
> st_intersects(site, reach, sparse=FALSE)
[,1]
[1,] FALSE
> st_intersects(site_snap, reach, sparse=FALSE)
[,1]
[1,] TRUE
Your original site is 10^-11 units from the line feature:
> st_distance(site$geom, reach$geom)
Units: m
[,1]
[1,] 9.550647e-11
Best Answer
Something like this: