shapefile - Expedite Near function in R for 10^9 coordinates?

Goal I want to find the geodesic distance between the closest hydrography feature to each of my lat/long coordinates.

I have hydrography data from New York state in a shapefile. I have many points, 10^9 number of coordinates in WGS1984 datum. There's about 20,000 features in my hydrography data which was originally in NAD1983.

I'm trying to find a way to expedite my code. This is what I've tried:

1) dist2Line from the package(geosphere): I keep getting the error: "returning Infno non-missing arguments to min". When I can get it to work (I haven't figured out why sometimes it works and sometimes it doesn't), it takes too long as well.

2) gDistance from package(rgeos): this also expects planar coordinates and takes too long (~43 sec for 300 pts).

3) st_nearest_feature from package (sf): This is fast (1 sec; 300 pts) but assumes the coordinates are planar. I get the index of the nearest hydro feature and I don't really know how to calculate distance easily from there (I suppose I could find a centroid for each hydro feature and calculate distance to that feature)?

4) st_nearest_points from package (sf): this takes too long (~100 sec, 300 pts) and assumes that coordinates are planar. It gives me linestring so I suppose I could calculate the distance between the two coordinates using a basic distance formula but the time & projection is not ideal.

I'm at a lost at how to proceed in a timely manner from here.

Potential Thoughts:

I've converted all of my hydrography features to WGS 1984 to work with it. Should I consider converting both files into planar coordinates? I'm not sure I understand well enough whether converting from unprojected to projected would cause any potential problems.

I just want the fastest way to process all of these points. I've been using fast_extract from package(velox) and (prioritizr) with raster data, and it's been incredibly fast. Something like that would be ideal.

I've tried to expedite this in ArcMap or R, and am open to either method. I've created a second question about ArcMap in the event somebody can think of a better me

Expedite Near function in ArcMap 10.7 for 10^9+ coordinates?