sf_st_distance_is_pretty_ok.R

library(sf)
options(stringsAsFactors = FALSE)

# 94720 points - 
# src: http://qldspatial.information.qld.gov.au/catalogue/custom/search.page?q=SALI_SITE
many_points <- st_read('D:/Spatial_data/soils/soil_survey_sites.shp') %>%
  st_transform(., crs = 3577) # always dist planar, kids

# divide in 2
many_points$nrow <- as.numeric(1:nrow(many_points))
points_odd <- many_points[many_points$nrow %% 2 == 0, ]
points_even <- many_points[setdiff(many_points$nrow, points_odd$nrow), ]

# demo
plot(st_geometry(points_odd[1:10, ]), pch = 19, col = 'blue')
plot(st_geometry(points_even[1:10, ]), add = TRUE, pch = 19, col = 'red')

t <- proc.time()
distmat <- st_distance(st_geometry(points_odd), 
                       st_geometry(points_even), 
                       by_element = TRUE)
t <- proc.time() - t
# ~30 seconds on an 8th gen Intel Core i7 with 8GB RAM
# returns a vector of distances in crs units (meters in this case); element[1]
# is the distance between point_odd[1,] and point_even[1,]