From 3f560365eb61546102f1322623cee3e6459e1fc3 Mon Sep 17 00:00:00 2001
From: olivroy <olivierroy71@hotmail.com>
Date: Mon, 2 Sep 2024 13:52:36 -0400
Subject: [PATCH] Fix #2428

---
 R/nearest.R              | 271 ++++++++++++++++++++-------------------
 man/st_nearest_points.Rd |  14 +-
 2 files changed, 143 insertions(+), 142 deletions(-)

diff --git a/R/nearest.R b/R/nearest.R
index 3a8728687..3cff5be78 100644
--- a/R/nearest.R
+++ b/R/nearest.R
@@ -1,135 +1,136 @@
-#' get nearest points between pairs of geometries
-#' 
-#' get nearest points between pairs of geometries
-#' @param x object of class \code{sfg}, \code{sfc} or \code{sf}
-#' @param y object of class \code{sfg}, \code{sfc} or \code{sf}
-#' @param pairwise logical; if \code{FALSE} (default) return nearest points between all pairs, if \code{TRUE}, return nearest points between subsequent pairs.
-#' @param ... ignored
-#' @seealso \link{st_nearest_feature} for finding the nearest feature
-#' @return an \link{sfc} object with all two-point \code{LINESTRING} geometries of point pairs from the first to the second geometry, of length x * y, with y cycling fastest. See examples for ideas how to convert these to \code{POINT} geometries.
-#' @details in case \code{x} lies inside \code{y}, when using S2, the end points 
-#' are on polygon boundaries, when using GEOS the end point are identical to \code{x}.
-#' @examples
-#' r = sqrt(2)/10
-#' pt1 = st_point(c(.1,.1))
-#' pt2 = st_point(c(.9,.9))
-#' pt3 = st_point(c(.9,.1))
-#' b1 = st_buffer(pt1, r)
-#' b2 = st_buffer(pt2, r)
-#' b3 = st_buffer(pt3, r)
-#' (ls0 = st_nearest_points(b1, b2)) # sfg
-#' (ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc
-#' plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green')
-#' plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue')
-#' plot(ls, add = TRUE, col = 'red')
-#' 
-#' nc = st_read(system.file("gpkg/nc.gpkg", package="sf"))
-#' plot(st_geometry(nc))
-#' ls = st_nearest_points(nc[1,], nc)
-#' plot(ls, col = 'red', add = TRUE)
-#' pts = st_cast(ls, "POINT") # gives all start & end points
-#' # starting, "from" points, corresponding to x:
-#' plot(pts[seq(1, 200, 2)], add = TRUE, col = 'blue')
-#' # ending, "to" points, corresponding to y:
-#' plot(pts[seq(2, 200, 2)], add = TRUE, col = 'green')
-#' 
-#' @export
-st_nearest_points = function(x, y, ...) UseMethod("st_nearest_points")
-
-#' @export
-#' @name st_nearest_points
-st_nearest_points.sfc = function(x, y, ..., pairwise = FALSE) {
-	stopifnot(st_crs(x) == st_crs(y))
-	longlat = isTRUE(st_is_longlat(x))
-	if (longlat && sf_use_s2()) {
-		ret = if (pairwise)
-				s2::s2_minimum_clearance_line_between(x, y)
-			else
-				do.call(c, lapply(x, s2::s2_minimum_clearance_line_between, y))
-		st_as_sfc(ret, crs = st_crs(x))
-	} else {
-		if (longlat)
-			message_longlat("st_nearest_points")
-		st_sfc(CPL_geos_nearest_points(x, st_geometry(y), pairwise), crs = st_crs(x))
-	}
-}
- 
-#' @export
-#' @name st_nearest_points
-st_nearest_points.sfg = function(x, y, ...) {
-	st_nearest_points(st_geometry(x), st_geometry(y), ...)
-}
-
-#' @export
-#' @name st_nearest_points
-st_nearest_points.sf = function(x, y, ...) {
-	st_nearest_points(st_geometry(x), st_geometry(y), ...)
-}
-
-#' get index of nearest feature
-#' 
-#' get index of nearest feature
-#' @param x object of class \code{sfg}, \code{sfc} or \code{sf}
-#' @param y object of class \code{sfg}, \code{sfc} or \code{sf}; if missing, features in \code{x} will be compared to all remaining features in \code{x}.
-#' @param ... ignored
-#' @param check_crs logical; should \code{x} and \code{y} be checked for CRS equality?
-#' @param longlat logical; does \code{x} have ellipsoidal coordinates?
-#' @return for each feature (geometry) in \code{x} the index of the nearest feature (geometry) in 
-#' set \code{y}, or in the remaining set of \code{x} if \code{y} is missing;
-#' empty geometries result in \code{NA} indexes
-#' @seealso \link{st_nearest_points} for finding the nearest points for pairs of feature geometries
-#' @export
-#' @examples
-#' ls1 = st_linestring(rbind(c(0,0), c(1,0)))
-#' ls2 = st_linestring(rbind(c(0,0.1), c(1,0.1)))
-#' ls3 = st_linestring(rbind(c(0,1), c(1,1)))
-#' (l = st_sfc(ls1, ls2, ls3))
-#' 
-#' p1 = st_point(c(0.1, -0.1))
-#' p2 = st_point(c(0.1, 0.11))
-#' p3 = st_point(c(0.1, 0.09))
-#' p4 = st_point(c(0.1, 0.9))
-#' 
-#' (p = st_sfc(p1, p2, p3, p4))
-#' try(st_nearest_feature(p, l))
-#' try(st_nearest_points(p, l[st_nearest_feature(p,l)], pairwise = TRUE))
-#' 
-#' r = sqrt(2)/10
-#' b1 = st_buffer(st_point(c(.1,.1)), r)
-#' b2 = st_buffer(st_point(c(.9,.9)), r)
-#' b3 = st_buffer(st_point(c(.9,.1)), r)
-#' circles = st_sfc(b1, b2, b3)
-#' plot(circles, col = NA, border = 2:4)
-#' pts = st_sfc(st_point(c(.3,.1)), st_point(c(.6,.2)), st_point(c(.6,.6)), st_point(c(.4,.8)))
-#' plot(pts, add = TRUE, col = 1)
-#' # draw points to nearest circle:
-#' nearest = try(st_nearest_feature(pts, circles))
-#' if (inherits(nearest, "try-error")) # GEOS 3.6.1 not available
-#'   nearest = c(1, 3, 2, 2)
-#' ls = st_nearest_points(pts, circles[nearest], pairwise = TRUE)
-#' plot(ls, col = 5:8, add = TRUE)
-#' # compute distance between pairs of nearest features:
-#' st_distance(pts, circles[nearest], by_element = TRUE)
-st_nearest_feature = function(x, y, ..., check_crs = TRUE, longlat = isTRUE(st_is_longlat(x))) {
-
-	if (missing(y)) { # https://github.com/r-spatial/s2/issues/111#issuecomment-835306261
-		longlat = force(longlat) # evaluate only once
-		x = st_geometry(x)
-		ind <- vapply(
-			seq_along(x),
-			function(i) st_nearest_feature(x[i], x[-i], check_crs = FALSE, longlat = longlat),
-			integer(1)
-		)
-		ifelse(ind >= seq_along(x), ind + 1, ind)
-	} else {
-		if (check_crs)
-			stopifnot(st_crs(x) == st_crs(y))
-		if (longlat && sf_use_s2())
-			s2::s2_closest_feature(x, y)
-		else {
-			if (longlat)
-				message_longlat("st_nearest_feature")
-			CPL_geos_nearest_feature(st_geometry(x), st_geometry(y))
-		}
-	}
-}
+#' Get nearest points between pairs of geometries
+#' 
+#' get nearest points between pairs of geometries
+#' @param x,y object of class \code{sfg}, \code{sfc} or \code{sf}
+#' @param pairwise logical; if \code{FALSE} (default) return nearest points between all pairs,
+#'   if \code{TRUE}, return nearest points between subsequent pairs.
+#' @param ... passed on to methods. Currently, only `pairwise` is implemented.
+#' @seealso \link{st_nearest_feature} for finding the nearest feature
+#' @return an \link{sfc} object with all two-point \code{LINESTRING} geometries of point pairs from the first to the second geometry, of length x * y, with y cycling fastest.
+#'   See examples for ideas how to convert these to \code{POINT} geometries.
+#' @details in case \code{x} lies inside \code{y}, when using S2, the end points 
+#' are on polygon boundaries, when using GEOS the end point are identical to \code{x}.
+#' @examples
+#' r = sqrt(2)/10
+#' pt1 = st_point(c(.1,.1))
+#' pt2 = st_point(c(.9,.9))
+#' pt3 = st_point(c(.9,.1))
+#' b1 = st_buffer(pt1, r)
+#' b2 = st_buffer(pt2, r)
+#' b3 = st_buffer(pt3, r)
+#' (ls0 = st_nearest_points(b1, b2)) # sfg
+#' (ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc
+#' plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green')
+#' plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue')
+#' plot(ls, add = TRUE, col = 'red')
+#' 
+#' nc = st_read(system.file("gpkg/nc.gpkg", package="sf"))
+#' plot(st_geometry(nc))
+#' ls = st_nearest_points(nc[1,], nc)
+#' plot(ls, col = 'red', add = TRUE)
+#' pts = st_cast(ls, "POINT") # gives all start & end points
+#' # starting, "from" points, corresponding to x:
+#' plot(pts[seq(1, 200, 2)], add = TRUE, col = 'blue')
+#' # ending, "to" points, corresponding to y:
+#' plot(pts[seq(2, 200, 2)], add = TRUE, col = 'green')
+#' 
+#' @export
+st_nearest_points = function(x, y, ...) UseMethod("st_nearest_points")
+
+#' @export
+#' @name st_nearest_points
+st_nearest_points.sfc = function(x, y, ..., pairwise = FALSE) {
+	stopifnot(st_crs(x) == st_crs(y))
+	longlat = isTRUE(st_is_longlat(x))
+	if (longlat && sf_use_s2()) {
+		ret = if (pairwise)
+				s2::s2_minimum_clearance_line_between(x, y)
+			else
+				do.call(c, lapply(x, s2::s2_minimum_clearance_line_between, y))
+		st_as_sfc(ret, crs = st_crs(x))
+	} else {
+		if (longlat)
+			message_longlat("st_nearest_points")
+		st_sfc(CPL_geos_nearest_points(x, st_geometry(y), pairwise), crs = st_crs(x))
+	}
+}
+ 
+#' @export
+#' @rdname st_nearest_points
+st_nearest_points.sfg = function(x, y, ...) {
+	st_nearest_points(st_geometry(x), st_geometry(y), ...)
+}
+
+#' @export
+#' @rdname st_nearest_points
+st_nearest_points.sf = function(x, y, ...) {
+	st_nearest_points(st_geometry(x), st_geometry(y), ...)
+}
+
+#' get index of nearest feature
+#' 
+#' get index of nearest feature
+#' @param x object of class \code{sfg}, \code{sfc} or \code{sf}
+#' @param y object of class \code{sfg}, \code{sfc} or \code{sf}; if missing, features in \code{x} will be compared to all remaining features in \code{x}.
+#' @param ... ignored
+#' @param check_crs logical; should \code{x} and \code{y} be checked for CRS equality?
+#' @param longlat logical; does \code{x} have ellipsoidal coordinates?
+#' @return for each feature (geometry) in \code{x} the index of the nearest feature (geometry) in 
+#' set \code{y}, or in the remaining set of \code{x} if \code{y} is missing;
+#' empty geometries result in \code{NA} indexes
+#' @seealso \link{st_nearest_points} for finding the nearest points for pairs of feature geometries
+#' @export
+#' @examples
+#' ls1 = st_linestring(rbind(c(0,0), c(1,0)))
+#' ls2 = st_linestring(rbind(c(0,0.1), c(1,0.1)))
+#' ls3 = st_linestring(rbind(c(0,1), c(1,1)))
+#' (l = st_sfc(ls1, ls2, ls3))
+#' 
+#' p1 = st_point(c(0.1, -0.1))
+#' p2 = st_point(c(0.1, 0.11))
+#' p3 = st_point(c(0.1, 0.09))
+#' p4 = st_point(c(0.1, 0.9))
+#' 
+#' (p = st_sfc(p1, p2, p3, p4))
+#' try(st_nearest_feature(p, l))
+#' try(st_nearest_points(p, l[st_nearest_feature(p,l)], pairwise = TRUE))
+#' 
+#' r = sqrt(2)/10
+#' b1 = st_buffer(st_point(c(.1,.1)), r)
+#' b2 = st_buffer(st_point(c(.9,.9)), r)
+#' b3 = st_buffer(st_point(c(.9,.1)), r)
+#' circles = st_sfc(b1, b2, b3)
+#' plot(circles, col = NA, border = 2:4)
+#' pts = st_sfc(st_point(c(.3,.1)), st_point(c(.6,.2)), st_point(c(.6,.6)), st_point(c(.4,.8)))
+#' plot(pts, add = TRUE, col = 1)
+#' # draw points to nearest circle:
+#' nearest = try(st_nearest_feature(pts, circles))
+#' if (inherits(nearest, "try-error")) # GEOS 3.6.1 not available
+#'   nearest = c(1, 3, 2, 2)
+#' ls = st_nearest_points(pts, circles[nearest], pairwise = TRUE)
+#' plot(ls, col = 5:8, add = TRUE)
+#' # compute distance between pairs of nearest features:
+#' st_distance(pts, circles[nearest], by_element = TRUE)
+st_nearest_feature = function(x, y, ..., check_crs = TRUE, longlat = isTRUE(st_is_longlat(x))) {
+
+	if (missing(y)) { # https://github.com/r-spatial/s2/issues/111#issuecomment-835306261
+		longlat = force(longlat) # evaluate only once
+		x = st_geometry(x)
+		ind <- vapply(
+			seq_along(x),
+			function(i) st_nearest_feature(x[i], x[-i], check_crs = FALSE, longlat = longlat),
+			integer(1)
+		)
+		ifelse(ind >= seq_along(x), ind + 1, ind)
+	} else {
+		if (check_crs)
+			stopifnot(st_crs(x) == st_crs(y))
+		if (longlat && sf_use_s2())
+			s2::s2_closest_feature(x, y)
+		else {
+			if (longlat)
+				message_longlat("st_nearest_feature")
+			CPL_geos_nearest_feature(st_geometry(x), st_geometry(y))
+		}
+	}
+}
diff --git a/man/st_nearest_points.Rd b/man/st_nearest_points.Rd
index 7fdfdeb1a..3e1f79802 100644
--- a/man/st_nearest_points.Rd
+++ b/man/st_nearest_points.Rd
@@ -5,7 +5,7 @@
 \alias{st_nearest_points.sfc}
 \alias{st_nearest_points.sfg}
 \alias{st_nearest_points.sf}
-\title{get nearest points between pairs of geometries}
+\title{Get nearest points between pairs of geometries}
 \usage{
 st_nearest_points(x, y, ...)
 
@@ -16,16 +16,16 @@ st_nearest_points(x, y, ...)
 \method{st_nearest_points}{sf}(x, y, ...)
 }
 \arguments{
-\item{x}{object of class \code{sfg}, \code{sfc} or \code{sf}}
+\item{x, y}{object of class \code{sfg}, \code{sfc} or \code{sf}}
 
-\item{y}{object of class \code{sfg}, \code{sfc} or \code{sf}}
+\item{...}{passed on to methods. Currently, only \code{pairwise} is implemented.}
 
-\item{...}{ignored}
-
-\item{pairwise}{logical; if \code{FALSE} (default) return nearest points between all pairs, if \code{TRUE}, return nearest points between subsequent pairs.}
+\item{pairwise}{logical; if \code{FALSE} (default) return nearest points between all pairs,
+if \code{TRUE}, return nearest points between subsequent pairs.}
 }
 \value{
-an \link{sfc} object with all two-point \code{LINESTRING} geometries of point pairs from the first to the second geometry, of length x * y, with y cycling fastest. See examples for ideas how to convert these to \code{POINT} geometries.
+an \link{sfc} object with all two-point \code{LINESTRING} geometries of point pairs from the first to the second geometry, of length x * y, with y cycling fastest.
+See examples for ideas how to convert these to \code{POINT} geometries.
 }
 \description{
 get nearest points between pairs of geometries