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?

  • You could give package geodist a try. It's supposed to be super fast and works with unprojected data IIRC. The caveat is that it only does point-to-point (I think). But maybe you could use all coordinate pairs of your line features as input and infer the feature by the line id (just a thought). – TimSalabim Dec 11 '19 at 17:30
  • 2
    @TimSalabim the nearest line vertex point might not belong to the nearest line segment. I think you have to compute the point-to-line distance and you can't fake it by only doing point to point. Densifying might help if you densify so all segments are smaller than the smallest real segment... but not convinced... – Spacedman Dec 11 '19 at 23:33
  • @Spacedman agreed. I don't think there is a straight-forward out-of-the-box solution for this. I guess an acceptable approach would be to first identify st_nearest_feature and then for each of the nearest features st_line_sample an appropriately high amount of points to do the actual distance calculation. As always, things depend on the context, hence no universal solution. – TimSalabim Dec 12 '19 at 9:44
  • Its never clear if any of these things is building a spatial index, or how long an overhead is spent building one. – Spacedman Dec 12 '19 at 9:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.