fixes #437

NEWS.md

@@ -1,5 +1,7 @@
 # version 0.5-3
 
+* `st_distance` gains a parameter `by_element` to obtain pairwise distances; #437
+
 * add the ability to `aggregate` using a simple feature `by` argument; #429
 
 * make the `op` argument to `[.sf` work

R/geom.R

@@ -120,14 +120,23 @@ st_geos_binop = function(op = "intersects", x, y, par = 0.0, pattern = NA_charac
 #' @param x object of class \code{sf}, \code{sfc} or \code{sfg}
 #' @param y object of class \code{sf}, \code{sfc} or \code{sfg}, defaults to \code{x}
 #' @param dist_fun function to be used for great circle distances of geographical coordinates; for unprojected (long/lat) data, this should be a distance function of package geosphere, or compatible to that; it defaults to \link[geosphere]{distGeo} in that case; for other data metric lengths are computed.
-#' @return st_distance returns a dense numeric matrix of dimension length(x) by length(y)
+#' @param by_element logical; if \code{TRUE}, return a vector with distance between the first elements of \code{x} and \code{y}, the second, etc. if \code{FALSE}, return the dense matrix with all pairwise distances.
+#' @return if \code{by_element} is \code{FALSE} a dense numeric matrix of dimension length(x) by length(y); otherwise a numeric vector of length \code{x} or \code{y}, the shorter one being recycled.
 #' @details function \code{dist_fun} should follow the pattern of the distance function \link[geosphere]{distGeo}: the first two arguments must be 2-column point matrices, the third the semi major axis (radius, in m), the third the ellipsoid flattening.
+#' @examples
+#' p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2)))
+#' st_distance(p, p)
+#' st_distance(p, p, by_element = TRUE)
 #' @export
-st_distance = function(x, y, dist_fun) {
+st_distance = function(x, y, dist_fun, by_element = FALSE) {
 	if (missing(y))
 		y = x
 	else
 		stopifnot(st_crs(x) == st_crs(y))
+
+	if (by_element)
+		return(mapply(st_distance, x, y, by_element = FALSE))
+
 	x = st_geometry(x)
 	y = st_geometry(y)
 	if (!is.na(st_crs(x)))

tests/geos.R

@@ -102,3 +102,7 @@ try(st_is_within_distance(nc_3857, nc_3857, 100000, sparse = FALSE))
 x = st_is_within_distance(nc_3857, nc_3857, 100000)
 y = st_is_within_distance(nc_3857, nc_3857, units::set_units(100, km))
 all.equal(x, y)
+
+p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2)))
+st_distance(p, p)
+st_distance(p, p, by_element = TRUE)

tests/geos.Rout.save

@@ -297,6 +297,15 @@ Error in CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern,  :
 > all.equal(x, y)
 [1] TRUE
 > 
+> p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2)))
+> st_distance(p, p)
+     [,1] [,2] [,3]
+[1,]    0    1    2
+[2,]    1    0    1
+[3,]    2    1    0
+> st_distance(p, p, by_element = TRUE)
+[1] 0 0 0
+> 
 > proc.time()
    user  system elapsed 
-  3.648   0.308   3.655 
+  3.764   0.244   3.756