Skip to content

Commit 3f56036

Browse files
olivroyolivroy
committed
Fix #2428 
1 parent ac32831 commit 3f56036

File tree

2 files changed

+143
-142
lines changed

2 files changed

+143
-142
lines changed

‎R/nearest.R

+136-135
Original file line numberDiff line numberDiff line change
@@ -1,135 +1,136 @@
1-
#' get nearest points between pairs of geometries
2-
#'
3-
#' get nearest points between pairs of geometries
4-
#' @param x object of class \code{sfg}, \code{sfc} or \code{sf}
5-
#' @param y object of class \code{sfg}, \code{sfc} or \code{sf}
6-
#' @param pairwise logical; if \code{FALSE} (default) return nearest points between all pairs, if \code{TRUE}, return nearest points between subsequent pairs.
7-
#' @param ... ignored
8-
#' @seealso \link{st_nearest_feature} for finding the nearest feature
9-
#' @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.
10-
#' @details in case \code{x} lies inside \code{y}, when using S2, the end points
11-
#' are on polygon boundaries, when using GEOS the end point are identical to \code{x}.
12-
#' @examples
13-
#' r = sqrt(2)/10
14-
#' pt1 = st_point(c(.1,.1))
15-
#' pt2 = st_point(c(.9,.9))
16-
#' pt3 = st_point(c(.9,.1))
17-
#' b1 = st_buffer(pt1, r)
18-
#' b2 = st_buffer(pt2, r)
19-
#' b3 = st_buffer(pt3, r)
20-
#' (ls0 = st_nearest_points(b1, b2)) # sfg
21-
#' (ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc
22-
#' plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green')
23-
#' plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue')
24-
#' plot(ls, add = TRUE, col = 'red')
25-
#'
26-
#' nc = st_read(system.file("gpkg/nc.gpkg", package="sf"))
27-
#' plot(st_geometry(nc))
28-
#' ls = st_nearest_points(nc[1,], nc)
29-
#' plot(ls, col = 'red', add = TRUE)
30-
#' pts = st_cast(ls, "POINT") # gives all start & end points
31-
#' # starting, "from" points, corresponding to x:
32-
#' plot(pts[seq(1, 200, 2)], add = TRUE, col = 'blue')
33-
#' # ending, "to" points, corresponding to y:
34-
#' plot(pts[seq(2, 200, 2)], add = TRUE, col = 'green')
35-
#'
36-
#' @export
37-
st_nearest_points = function(x, y, ...) UseMethod("st_nearest_points")
38-
39-
#' @export
40-
#' @name st_nearest_points
41-
st_nearest_points.sfc = function(x, y, ..., pairwise = FALSE) {
42-
stopifnot(st_crs(x) == st_crs(y))
43-
longlat = isTRUE(st_is_longlat(x))
44-
if (longlat && sf_use_s2()) {
45-
ret = if (pairwise)
46-
s2::s2_minimum_clearance_line_between(x, y)
47-
else
48-
do.call(c, lapply(x, s2::s2_minimum_clearance_line_between, y))
49-
st_as_sfc(ret, crs = st_crs(x))
50-
} else {
51-
if (longlat)
52-
message_longlat("st_nearest_points")
53-
st_sfc(CPL_geos_nearest_points(x, st_geometry(y), pairwise), crs = st_crs(x))
54-
}
55-
}
56-
57-
#' @export
58-
#' @name st_nearest_points
59-
st_nearest_points.sfg = function(x, y, ...) {
60-
st_nearest_points(st_geometry(x), st_geometry(y), ...)
61-
}
62-
63-
#' @export
64-
#' @name st_nearest_points
65-
st_nearest_points.sf = function(x, y, ...) {
66-
st_nearest_points(st_geometry(x), st_geometry(y), ...)
67-
}
68-
69-
#' get index of nearest feature
70-
#'
71-
#' get index of nearest feature
72-
#' @param x object of class \code{sfg}, \code{sfc} or \code{sf}
73-
#' @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}.
74-
#' @param ... ignored
75-
#' @param check_crs logical; should \code{x} and \code{y} be checked for CRS equality?
76-
#' @param longlat logical; does \code{x} have ellipsoidal coordinates?
77-
#' @return for each feature (geometry) in \code{x} the index of the nearest feature (geometry) in
78-
#' set \code{y}, or in the remaining set of \code{x} if \code{y} is missing;
79-
#' empty geometries result in \code{NA} indexes
80-
#' @seealso \link{st_nearest_points} for finding the nearest points for pairs of feature geometries
81-
#' @export
82-
#' @examples
83-
#' ls1 = st_linestring(rbind(c(0,0), c(1,0)))
84-
#' ls2 = st_linestring(rbind(c(0,0.1), c(1,0.1)))
85-
#' ls3 = st_linestring(rbind(c(0,1), c(1,1)))
86-
#' (l = st_sfc(ls1, ls2, ls3))
87-
#'
88-
#' p1 = st_point(c(0.1, -0.1))
89-
#' p2 = st_point(c(0.1, 0.11))
90-
#' p3 = st_point(c(0.1, 0.09))
91-
#' p4 = st_point(c(0.1, 0.9))
92-
#'
93-
#' (p = st_sfc(p1, p2, p3, p4))
94-
#' try(st_nearest_feature(p, l))
95-
#' try(st_nearest_points(p, l[st_nearest_feature(p,l)], pairwise = TRUE))
96-
#'
97-
#' r = sqrt(2)/10
98-
#' b1 = st_buffer(st_point(c(.1,.1)), r)
99-
#' b2 = st_buffer(st_point(c(.9,.9)), r)
100-
#' b3 = st_buffer(st_point(c(.9,.1)), r)
101-
#' circles = st_sfc(b1, b2, b3)
102-
#' plot(circles, col = NA, border = 2:4)
103-
#' pts = st_sfc(st_point(c(.3,.1)), st_point(c(.6,.2)), st_point(c(.6,.6)), st_point(c(.4,.8)))
104-
#' plot(pts, add = TRUE, col = 1)
105-
#' # draw points to nearest circle:
106-
#' nearest = try(st_nearest_feature(pts, circles))
107-
#' if (inherits(nearest, "try-error")) # GEOS 3.6.1 not available
108-
#' nearest = c(1, 3, 2, 2)
109-
#' ls = st_nearest_points(pts, circles[nearest], pairwise = TRUE)
110-
#' plot(ls, col = 5:8, add = TRUE)
111-
#' # compute distance between pairs of nearest features:
112-
#' st_distance(pts, circles[nearest], by_element = TRUE)
113-
st_nearest_feature = function(x, y, ..., check_crs = TRUE, longlat = isTRUE(st_is_longlat(x))) {
114-
115-
if (missing(y)) { # https://github.com/r-spatial/s2/issues/111#issuecomment-835306261
116-
longlat = force(longlat) # evaluate only once
117-
x = st_geometry(x)
118-
ind <- vapply(
119-
seq_along(x),
120-
function(i) st_nearest_feature(x[i], x[-i], check_crs = FALSE, longlat = longlat),
121-
integer(1)
122-
)
123-
ifelse(ind >= seq_along(x), ind + 1, ind)
124-
} else {
125-
if (check_crs)
126-
stopifnot(st_crs(x) == st_crs(y))
127-
if (longlat && sf_use_s2())
128-
s2::s2_closest_feature(x, y)
129-
else {
130-
if (longlat)
131-
message_longlat("st_nearest_feature")
132-
CPL_geos_nearest_feature(st_geometry(x), st_geometry(y))
133-
}
134-
}
135-
}
1+
#' Get nearest points between pairs of geometries
2+
#'
3+
#' get nearest points between pairs of geometries
4+
#' @param x,y object of class \code{sfg}, \code{sfc} or \code{sf}
5+
#' @param pairwise logical; if \code{FALSE} (default) return nearest points between all pairs,
6+
#' if \code{TRUE}, return nearest points between subsequent pairs.
7+
#' @param ... passed on to methods. Currently, only `pairwise` is implemented.
8+
#' @seealso \link{st_nearest_feature} for finding the nearest feature
9+
#' @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.
10+
#' See examples for ideas how to convert these to \code{POINT} geometries.
11+
#' @details in case \code{x} lies inside \code{y}, when using S2, the end points
12+
#' are on polygon boundaries, when using GEOS the end point are identical to \code{x}.
13+
#' @examples
14+
#' r = sqrt(2)/10
15+
#' pt1 = st_point(c(.1,.1))
16+
#' pt2 = st_point(c(.9,.9))
17+
#' pt3 = st_point(c(.9,.1))
18+
#' b1 = st_buffer(pt1, r)
19+
#' b2 = st_buffer(pt2, r)
20+
#' b3 = st_buffer(pt3, r)
21+
#' (ls0 = st_nearest_points(b1, b2)) # sfg
22+
#' (ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc
23+
#' plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green')
24+
#' plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue')
25+
#' plot(ls, add = TRUE, col = 'red')
26+
#'
27+
#' nc = st_read(system.file("gpkg/nc.gpkg", package="sf"))
28+
#' plot(st_geometry(nc))
29+
#' ls = st_nearest_points(nc[1,], nc)
30+
#' plot(ls, col = 'red', add = TRUE)
31+
#' pts = st_cast(ls, "POINT") # gives all start & end points
32+
#' # starting, "from" points, corresponding to x:
33+
#' plot(pts[seq(1, 200, 2)], add = TRUE, col = 'blue')
34+
#' # ending, "to" points, corresponding to y:
35+
#' plot(pts[seq(2, 200, 2)], add = TRUE, col = 'green')
36+
#'
37+
#' @export
38+
st_nearest_points = function(x, y, ...) UseMethod("st_nearest_points")
39+
40+
#' @export
41+
#' @name st_nearest_points
42+
st_nearest_points.sfc = function(x, y, ..., pairwise = FALSE) {
43+
stopifnot(st_crs(x) == st_crs(y))
44+
longlat = isTRUE(st_is_longlat(x))
45+
if (longlat && sf_use_s2()) {
46+
ret = if (pairwise)
47+
s2::s2_minimum_clearance_line_between(x, y)
48+
else
49+
do.call(c, lapply(x, s2::s2_minimum_clearance_line_between, y))
50+
st_as_sfc(ret, crs = st_crs(x))
51+
} else {
52+
if (longlat)
53+
message_longlat("st_nearest_points")
54+
st_sfc(CPL_geos_nearest_points(x, st_geometry(y), pairwise), crs = st_crs(x))
55+
}
56+
}
57+
58+
#' @export
59+
#' @rdname st_nearest_points
60+
st_nearest_points.sfg = function(x, y, ...) {
61+
st_nearest_points(st_geometry(x), st_geometry(y), ...)
62+
}
63+
64+
#' @export
65+
#' @rdname st_nearest_points
66+
st_nearest_points.sf = function(x, y, ...) {
67+
st_nearest_points(st_geometry(x), st_geometry(y), ...)
68+
}
69+
70+
#' get index of nearest feature
71+
#'
72+
#' get index of nearest feature
73+
#' @param x object of class \code{sfg}, \code{sfc} or \code{sf}
74+
#' @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}.
75+
#' @param ... ignored
76+
#' @param check_crs logical; should \code{x} and \code{y} be checked for CRS equality?
77+
#' @param longlat logical; does \code{x} have ellipsoidal coordinates?
78+
#' @return for each feature (geometry) in \code{x} the index of the nearest feature (geometry) in
79+
#' set \code{y}, or in the remaining set of \code{x} if \code{y} is missing;
80+
#' empty geometries result in \code{NA} indexes
81+
#' @seealso \link{st_nearest_points} for finding the nearest points for pairs of feature geometries
82+
#' @export
83+
#' @examples
84+
#' ls1 = st_linestring(rbind(c(0,0), c(1,0)))
85+
#' ls2 = st_linestring(rbind(c(0,0.1), c(1,0.1)))
86+
#' ls3 = st_linestring(rbind(c(0,1), c(1,1)))
87+
#' (l = st_sfc(ls1, ls2, ls3))
88+
#'
89+
#' p1 = st_point(c(0.1, -0.1))
90+
#' p2 = st_point(c(0.1, 0.11))
91+
#' p3 = st_point(c(0.1, 0.09))
92+
#' p4 = st_point(c(0.1, 0.9))
93+
#'
94+
#' (p = st_sfc(p1, p2, p3, p4))
95+
#' try(st_nearest_feature(p, l))
96+
#' try(st_nearest_points(p, l[st_nearest_feature(p,l)], pairwise = TRUE))
97+
#'
98+
#' r = sqrt(2)/10
99+
#' b1 = st_buffer(st_point(c(.1,.1)), r)
100+
#' b2 = st_buffer(st_point(c(.9,.9)), r)
101+
#' b3 = st_buffer(st_point(c(.9,.1)), r)
102+
#' circles = st_sfc(b1, b2, b3)
103+
#' plot(circles, col = NA, border = 2:4)
104+
#' pts = st_sfc(st_point(c(.3,.1)), st_point(c(.6,.2)), st_point(c(.6,.6)), st_point(c(.4,.8)))
105+
#' plot(pts, add = TRUE, col = 1)
106+
#' # draw points to nearest circle:
107+
#' nearest = try(st_nearest_feature(pts, circles))
108+
#' if (inherits(nearest, "try-error")) # GEOS 3.6.1 not available
109+
#' nearest = c(1, 3, 2, 2)
110+
#' ls = st_nearest_points(pts, circles[nearest], pairwise = TRUE)
111+
#' plot(ls, col = 5:8, add = TRUE)
112+
#' # compute distance between pairs of nearest features:
113+
#' st_distance(pts, circles[nearest], by_element = TRUE)
114+
st_nearest_feature = function(x, y, ..., check_crs = TRUE, longlat = isTRUE(st_is_longlat(x))) {
115+
116+
if (missing(y)) { # https://github.com/r-spatial/s2/issues/111#issuecomment-835306261
117+
longlat = force(longlat) # evaluate only once
118+
x = st_geometry(x)
119+
ind <- vapply(
120+
seq_along(x),
121+
function(i) st_nearest_feature(x[i], x[-i], check_crs = FALSE, longlat = longlat),
122+
integer(1)
123+
)
124+
ifelse(ind >= seq_along(x), ind + 1, ind)
125+
} else {
126+
if (check_crs)
127+
stopifnot(st_crs(x) == st_crs(y))
128+
if (longlat && sf_use_s2())
129+
s2::s2_closest_feature(x, y)
130+
else {
131+
if (longlat)
132+
message_longlat("st_nearest_feature")
133+
CPL_geos_nearest_feature(st_geometry(x), st_geometry(y))
134+
}
135+
}
136+
}

‎man/st_nearest_points.Rd

+7-7
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)