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