1#' @name geos_query
2#' @export
3#' @return st_is_simple returns a logical vector, indicating for each geometry whether it is simple (e.g., not self-intersecting)
4#' @examples
5#' ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1)))
6#' st_is_simple(st_sfc(ls, st_point(c(0,0))))
7st_is_simple = function(x) CPL_geos_is_simple(st_geometry(x))
8
9#' @name geos_query
10#' @export
11#' @return st_is_empty returns for each geometry whether it is empty
12#' @examples
13#' ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1)))
14#' st_is_empty(st_sfc(ls, st_point(), st_linestring()))
15st_is_empty = function(x) CPL_geos_is_empty(st_geometry(x))
16
17is_symmetric = function(operation, pattern) {
18	if (!is.na(pattern)) {
19		m = matrix(sapply(1:9, function(i) substr(pattern, i, i)), 3, 3)
20		isTRUE(all(m == t(m)))
21	} else
22		isTRUE(operation %in% c("intersects", "touches", "overlaps", "disjoint", "equals"))
23}
24
25# binary, interfaced through GEOS or S2:
26
27# [1] X "s2_contains_matrix"    X "s2_covered_by_matrix"
28# [3] X "s2_covers_matrix"      X "s2_disjoint_matrix"
29# [5] X "s2_distance_matrix"    X "s2_dwithin_matrix"
30# [7] X "s2_equals_matrix"      X "s2_intersects_matrix"
31# [9] "s2_max_distance_matrix"    "s2_may_intersect_matrix"
32#[11] X "s2_touches_matrix"     X "s2_within_matrix"
33
34# returning matrix, distance or relation string -- the work horse is:
35st_geos_binop = function(op, x, y, par = 0.0, pattern = NA_character_,
36		sparse = TRUE, prepared = FALSE, model = "closed", ...) {
37	if (missing(y))
38		y = x
39	else if (inherits(x, c("sf", "sfc")) && inherits(y, c("sf", "sfc")))
40		stopifnot(st_crs(x) == st_crs(y))
41	longlat = inherits(x, "s2geography") || isTRUE(st_is_longlat(x))
42	if (longlat && sf_use_s2() && op %in% c("intersects", "contains", "within",
43			"covers", "covered_by", "disjoint", "equals", "touches")) {
44		fn = get(paste0("s2_", op, "_matrix"), envir = getNamespace("s2")) # get op function
45		lst = fn(x, y, s2::s2_options(model = model, ...)) # call function
46		id = if (is.null(row.names(x)))
47				as.character(seq_along(lst))
48			else
49				row.names(x)
50		sgbp(lst, predicate = op, region.id = id, ncol = length(st_geometry(y)), sparse)
51	} else {
52		if (longlat && !(op %in% c("equals", "equals_exact")))
53			message_longlat(paste0("st_", op))
54		if (prepared && is_symmetric(op, pattern) &&
55				length(dx <- st_dimension(x)) && length(dy <- st_dimension(y)) &&
56				isTRUE(all(dx == 0)) && isTRUE(all(dy == 2))) {
57			t(st_geos_binop(op, y, x, par = par, pattern = pattern, sparse = sparse,
58				prepared = prepared))
59		} else {
60			ret = CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern, prepared)
61			if (length(ret) == 0 || is.null(dim(ret[[1]]))) {
62				id = if (is.null(row.names(x)))
63						as.character(seq_along(ret))
64					else
65						row.names(x)
66				sgbp(ret, predicate = op, region.id = id, ncol = length(st_geometry(y)), sparse)
67			} else # CPL_geos_binop returned a matrix, e.g. from op = "relate"
68				ret[[1]]
69		}
70	}
71}
72
73#' Compute DE9-IM relation between pairs of geometries, or match it to a given pattern
74#'
75#' Compute DE9-IM relation between pairs of geometries, or match it to a given pattern
76#' @param x object of class \code{sf}, \code{sfc} or \code{sfg}
77#' @param y object of class \code{sf}, \code{sfc} or \code{sfg}
78#' @param pattern character; define the pattern to match to, see details.
79#' @param sparse logical; should a sparse matrix be returned (TRUE) or a dense matrix?
80#' @return In case \code{pattern} is not given, \code{st_relate} returns a dense \code{character} matrix; element [i,j] has nine characters, referring to the DE9-IM relationship between x[i] and y[j], encoded as IxIy,IxBy,IxEy,BxIy,BxBy,BxEy,ExIy,ExBy,ExEy where I refers to interior, B to boundary, and E to exterior, and e.g. BxIy the dimensionality of the intersection of the the boundary of x[i] and the interior of y[j], which is one of {0,1,2,F}, digits denoting dimensionality, F denoting not intersecting. When \code{pattern} is given, a dense logical matrix or sparse index list returned with matches to the given pattern; see \link{st_intersection} for a description of the returned matrix or list. See also \url{https://en.wikipedia.org/wiki/DE-9IM} for further explanation.
81#' @export
82#' @examples
83#' p1 = st_point(c(0,0))
84#' p2 = st_point(c(2,2))
85#' pol1 = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0)))) - 0.5
86#' pol2 = pol1 + 1
87#' pol3 = pol1 + 2
88#' st_relate(st_sfc(p1, p2), st_sfc(pol1, pol2, pol3))
89#' sfc = st_sfc(st_point(c(0,0)), st_point(c(3,3)))
90#' grd = st_make_grid(sfc, n = c(3,3))
91#' st_intersects(grd)
92#' st_relate(grd, pattern = "****1****") # sides, not corners, internals
93#' st_relate(grd, pattern = "****0****") # only corners touch
94#' st_rook = function(a, b = a) st_relate(a, b, pattern = "F***1****")
95#' st_rook(grd)
96#' # queen neighbours, see \url{https://github.com/r-spatial/sf/issues/234#issuecomment-300511129}
97#' st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****")
98st_relate	= function(x, y, pattern = NA_character_, sparse = !is.na(pattern)) {
99	if (!is.na(pattern)) {
100		stopifnot(is.character(pattern) && length(pattern) == 1 && nchar(pattern) == 9)
101		st_geos_binop("relate_pattern", x, y, pattern = pattern, sparse = sparse)
102	} else
103		st_geos_binop("relate", x, y, sparse = FALSE)
104}
105
106#' Geometric binary predicates on pairs of simple feature geometry sets
107#'
108#' Geometric binary predicates on pairs of simple feature geometry sets
109#' @name geos_binary_pred
110#' @param x object of class \code{sf}, \code{sfc} or \code{sfg}
111#' @param y object of class \code{sf}, \code{sfc} or \code{sfg}; if missing, \code{x} is used
112#' @param sparse logical; should a sparse index list be returned (TRUE) or a dense logical matrix? See below.
113#' @param ... passed on to \link[s2]{s2_options}
114#' @param prepared logical; prepare geometry for x, before looping over y? See Details.
115#' @details If \code{prepared} is \code{TRUE}, and \code{x} contains POINT geometries and \code{y} contains polygons, then the polygon geometries are prepared, rather than the points.
116#' @return If \code{sparse=FALSE}, \code{st_predicate} (with \code{predicate} e.g. "intersects") returns a dense logical matrix with element \code{i,j} \code{TRUE} when \code{predicate(x[i], y[j])} (e.g., when geometry of feature i and j intersect); if \code{sparse=TRUE}, an object of class \code{\link{sgbp}} with a sparse list representation of the same matrix, with list element \code{i} an integer vector with all indices j for which \code{predicate(x[i],y[j])} is \code{TRUE} (and hence a zero-length integer vector if none of them is \code{TRUE}). From the dense matrix, one can find out if one or more elements intersect by \code{apply(mat, 1, any)}, and from the sparse list by \code{lengths(lst) > 0}, see examples below.
117#' @details For most predicates, a spatial index is built on argument \code{x}; see \url{https://r-spatial.org/r/2017/06/22/spatial-index.html}.
118#' Specifically, \code{st_intersects}, \code{st_disjoint}, \code{st_touches} \code{st_crosses}, \code{st_within}, \code{st_contains}, \code{st_contains_properly}, \code{st_overlaps}, \code{st_equals}, \code{st_covers} and \code{st_covered_by} all build spatial indexes for more efficient geometry calculations. \code{st_relate}, \code{st_equals_exact}, and do not; \code{st_is_within_distance} uses a spatial index for geographic coordinates when \code{sf_use_s2()} is true.
119#'
120#' If \code{y} is missing, `st_predicate(x, x)` is effectively called, and a square matrix is returned with diagonal elements `st_predicate(x[i], x[i])`.
121#'
122#' Sparse geometry binary predicate (\code{\link{sgbp}}) lists have the following attributes: \code{region.id} with the \code{row.names} of \code{x} (if any, else \code{1:n}), \code{ncol} with the number of features in \code{y}, and \code{predicate} with the name of the predicate used.
123#'
124#' @note For intersection on pairs of simple feature geometries, use
125#' the function \code{\link{st_intersection}} instead of \code{st_intersects}.
126#'
127#' @examples
128#' pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 2.5)))
129#' pol = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0))))
130#' (lst = st_intersects(pts, pol))
131#' (mat = st_intersects(pts, pol, sparse = FALSE))
132#' # which points fall inside a polygon?
133#' apply(mat, 1, any)
134#' lengths(lst) > 0
135#' # which points fall inside the first polygon?
136#' st_intersects(pol, pts)[[1]]
137#' @export
138st_intersects	= function(x, y, sparse = TRUE, ...) UseMethod("st_intersects")
139
140#' @export
141st_intersects.sfc = function(x, y, sparse = TRUE, prepared = TRUE, ...)
142	st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, ...)
143
144#' @export
145st_intersects.sf = function(x, y, sparse = TRUE, prepared = TRUE, ...)
146	st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, ...)
147
148#' @export
149st_intersects.sfg = function(x, y, sparse = TRUE, prepared = TRUE, ...)
150	st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, ...)
151
152
153#' @name geos_binary_pred
154#' @export
155st_disjoint		= function(x, y = x, sparse = TRUE, prepared = TRUE) {
156	# st_geos_binop("disjoint", x, y, sparse = sparse, prepared = prepared) -> didn't use STRtree
157	int = st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared)
158	# disjoint = !intersects :
159	if (sparse)
160		sgbp(lapply(int, function(g) setdiff(seq_along(st_geometry(y)), g)),
161			predicate = "disjoint",
162			ncol = attr(int, "ncol"),
163			region.id = attr(int, "region.id"))
164	else
165		!int
166}
167
168#' @name geos_binary_pred
169#' @export
170st_touches		= function(x, y, sparse = TRUE, prepared = TRUE, ...)
171	st_geos_binop("touches", x, y, sparse = sparse, prepared = prepared, ...)
172
173#' @name geos_binary_pred
174#' @export
175st_crosses		= function(x, y, sparse = TRUE, prepared = TRUE, ...)
176	st_geos_binop("crosses", x, y, sparse = sparse, prepared = prepared, ...)
177
178#' @name geos_binary_pred
179#' @export
180st_within	= function(x, y, sparse = TRUE, prepared = TRUE, ...)
181	st_geos_binop("within", x, y, sparse = sparse, prepared = prepared, ...)
182
183#' @name geos_binary_pred
184#' @param model character; polygon/polyline model; one of
185#' "open", "semi-open" or "closed"; see Details.
186#' @details for \code{model}, see https://github.com/r-spatial/s2/issues/32
187#' @export
188st_contains	= function(x, y, sparse = TRUE, prepared = TRUE, ..., model = "open")
189	st_geos_binop("contains", x, y, sparse = sparse, prepared = prepared, ..., model = model)
190
191#' @name geos_binary_pred
192#' @export
193#' @details `st_contains_properly(A,B)` is true if A intersects B's interior, but not its edges or exterior; A contains A, but A does not properly contain A.
194#'
195#' See also \link{st_relate} and \url{https://en.wikipedia.org/wiki/DE-9IM} for a more detailed description of the underlying algorithms.
196st_contains_properly = function(x, y, sparse = TRUE, prepared = TRUE, ...) {
197	if (! prepared)
198		stop("non-prepared geometries not supported for st_contains_properly")
199	st_geos_binop("contains_properly", x, y, sparse = sparse, prepared = TRUE, ...)
200}
201
202#' @name geos_binary_pred
203#' @export
204st_overlaps		= function(x, y, sparse = TRUE, prepared = TRUE, ...)
205	st_geos_binop("overlaps", x, y, sparse = sparse, prepared = prepared, ...)
206
207#' @name geos_binary_pred
208#' @export
209st_equals		= function(x, y, sparse = TRUE, prepared = FALSE, ...) {
210	if (prepared)
211		stop("prepared geometries not supported for st_equals")
212	st_geos_binop("equals", x, y, sparse = sparse, ...)
213}
214
215#' @name geos_binary_pred
216#' @export
217st_covers = function(x, y, sparse = TRUE, prepared = TRUE, ..., model = "closed")
218	st_geos_binop("covers", x, y, sparse = sparse, prepared = prepared, ..., model = model)
219
220
221#' @name geos_binary_pred
222#' @export
223st_covered_by = function(x, y = x, sparse = TRUE, prepared = TRUE, ..., model = "closed")
224	st_geos_binop("covered_by", x, y, sparse = sparse, prepared = prepared, ...)
225
226
227#' @name geos_binary_pred
228#' @export
229#' @param par numeric; parameter used for "equals_exact" (margin);
230#' @details \code{st_equals_exact} returns true for two geometries of the same type and their vertices corresponding by index are equal up to a specified tolerance.
231st_equals_exact = function(x, y, par, sparse = TRUE, prepared = FALSE, ...) {
232	if (prepared)
233		stop("prepared geometries not supported for st_equals_exact")
234	st_geos_binop("equals_exact", x, y, par = par, sparse = sparse, ...)
235}
236
237#' @name geos_binary_pred
238#' @export
239#' @param dist distance threshold; geometry indexes with distances smaller or equal to this value are returned; numeric value or units value having distance units.
240st_is_within_distance = function(x, y = x, dist, sparse = TRUE, ...) {
241
242	ret = if (isTRUE(st_is_longlat(x))) {
243			units(dist) = as_units("m") # might convert
244			r = if (sf_use_s2()) {
245				if (inherits(dist, "units"))
246					dist = drop_units(dist)
247				s2::s2_dwithin_matrix(x, y, dist, ...)
248			} else {
249				if (!requireNamespace("lwgeom", quietly = TRUE) ||
250						utils::packageVersion("lwgeom") <= "0.1-2")
251					stop("lwgeom > 0.1-2 required: install first?")
252				lwgeom::st_geod_distance(x, y, tolerance = dist, sparse = TRUE)
253			}
254			sgbp(r, predicate = "is_within_distance", region.id = seq_along(x),
255				ncol = length(st_geometry(y)))
256		} else {
257			if (! is.na(st_crs(x)))
258				units(dist) = crs_parameters(st_crs(x))$ud_unit # might convert
259			st_geos_binop("is_within_distance", x, y, par = dist, sparse = sparse, ...)
260		}
261	if (!sparse)
262		as.matrix(ret)
263	else
264		ret
265}
266