1#' @export 2#' @name st_sample 3st_sample = function(x, size, ...) UseMethod("st_sample") 4 5#' sample points on or in (sets of) spatial features 6#' 7#' Sample points on or in (sets of) spatial features. 8#' By default, returns a pre-specified number of points that is equal to 9#' \code{size} (if \code{type = "random"} and \code{exact = TRUE}) or an approximation of 10#' \code{size} otherwise. \code{spatstat} methods are 11#' interfaced and do not use the \code{size} argument, see examples. 12#' 13#' The function is vectorised: it samples \code{size} points across all geometries in 14#' the object if \code{size} is a single number, or the specified number of points 15#' in each feature if \code{size} is a vector of integers equal in length to the geometry 16#' of \code{x}. 17#' 18#' @param x object of class \code{sf} or \code{sfc} 19#' @param size sample size(s) requested; either total size, or a numeric vector with sample sizes for each feature geometry. When sampling polygons, the returned sampling size may differ from the requested size, as the bounding box is sampled, and sampled points intersecting the polygon are returned. 20#' @param warn_if_not_integer logical; if \code{FALSE} then no warning is emitted if \code{size} is not an integer 21#' @param ... passed on to \link[base]{sample} for \code{multipoint} sampling, or to \code{spatstat} functions for spatstat sampling types (see details) 22#' @param type character; indicates the spatial sampling type; one of \code{random}, \code{hexagonal} (triangular really), \code{regular}, 23#' or one of the \code{spatstat} methods such as \code{Thomas} for calling \code{spatstat.core::rThomas} (see Details). 24#' @param exact logical; should the length of output be exactly 25#' @param by_polygon logical; for \code{MULTIPOLYGON} geometries, should the effort be split by \code{POLYGON}? See https://github.com/r-spatial/sf/issues/1480 26#' the same as specified by \code{size}? \code{TRUE} by default. Only applies to polygons, and 27#' when \code{type = "random"}. 28#' @return an \code{sfc} object containing the sampled \code{POINT} geometries 29#' @details if \code{x} has dimension 2 (polygons) and geographical coordinates (long/lat), uniform random sampling on the sphere is applied, see e.g. \url{http://mathworld.wolfram.com/SpherePointPicking.html} 30#' 31#' For \code{regular} or \code{hexagonal} sampling of polygons, the resulting size is only an approximation. 32#' 33#' As parameter called \code{offset} can be passed to control ("fix") regular or hexagonal sampling: for polygons a length 2 numeric vector (by default: a random point from \code{st_bbox(x)}); for lines use a number like \code{runif(1)}. 34#' 35#' Sampling methods from package \code{spatstat} are interfaced (see examples), and need their own parameters to be set. 36#' For instance, to use \code{spatstat.core::rThomas()}, set \code{type = "Thomas"}. 37#' @examples 38#' nc = st_read(system.file("shape/nc.shp", package="sf")) 39#' p1 = st_sample(nc[1:3, ], 6) 40#' p2 = st_sample(nc[1:3, ], 1:3) 41#' plot(st_geometry(nc)[1:3]) 42#' plot(p1, add = TRUE) 43#' plot(p2, add = TRUE, pch = 2) 44#' x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,90),c(0,90),c(0,0)))), crs = st_crs(4326)) 45#' plot(x, axes = TRUE, graticule = TRUE) 46#' if (sf_extSoftVersion()["proj.4"] >= "4.9.0") 47#' plot(p <- st_sample(x, 1000), add = TRUE) 48#' x2 = st_transform(st_segmentize(x, 1e4), st_crs("+proj=ortho +lat_0=30 +lon_0=45")) 49#' g = st_transform(st_graticule(), st_crs("+proj=ortho +lat_0=30 +lon_0=45")) 50#' plot(x2, graticule = g) 51#' if (sf_extSoftVersion()["proj.4"] >= "4.9.0") { 52#' p2 = st_transform(p, st_crs("+proj=ortho +lat_0=30 +lon_0=45")) 53#' plot(p2, add = TRUE) 54#' } 55#' x = st_sfc(st_polygon(list(rbind(c(0,0),c(90,0),c(90,10),c(0,90),c(0,0))))) # NOT long/lat: 56#' plot(x) 57#' p_exact = st_sample(x, 1000, exact = TRUE) 58#' p_not_exact = st_sample(x, 1000, exact = FALSE) 59#' length(p_exact); length(p_not_exact) 60#' plot(st_sample(x, 1000), add = TRUE) 61#' x = st_sfc(st_polygon(list(rbind(c(-180,-90),c(180,-90),c(180,90),c(-180,90),c(-180,-90)))), 62#' crs=st_crs(4326)) 63#' # FIXME: 64#' #if (sf_extSoftVersion()["proj.4"] >= "4.9.0") { 65#' # p = st_sample(x, 1000) 66#' # st_sample(p, 3) 67#' #} 68#' # hexagonal: 69#' sfc = st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))) 70#' plot(sfc) 71#' h = st_sample(sfc, 100, type = "hexagonal") 72#' h1 = st_sample(sfc, 100, type = "hexagonal") 73#' plot(h, add = TRUE) 74#' plot(h1, col = 'red', add = TRUE) 75#' c(length(h), length(h1)) # approximate! 76#' pt = st_multipoint(matrix(1:20,,2)) 77#' ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), 78#' st_linestring(rbind(c(0,0),c(.1,0))), 79#' st_linestring(rbind(c(0,1),c(.1,1))), 80#' st_linestring(rbind(c(2,2),c(2,2.00001)))) 81#' st_sample(ls, 80) 82#' plot(st_sample(ls, 80)) 83#' # spatstat example: 84#' if (require(spatstat.core)) { 85#' x <- sf::st_sfc(sf::st_polygon(list(rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 0))))) 86#' # for spatstat.core::rThomas(), set type = "Thomas": 87#' pts <- st_sample(x, kappa = 1, mu = 10, scale = 0.1, type = "Thomas") 88#' } 89#' @export 90#' @name st_sample 91st_sample.sf = function(x, size, ...) st_sample(st_geometry(x), size, ...) 92 93#' @export 94#' @name st_sample 95st_sample.sfc = function(x, size, ..., type = "random", exact = TRUE, warn_if_not_integer = TRUE, 96 by_polygon = FALSE) { 97 98 if (!missing(size) && warn_if_not_integer && any(size %% 1 != 0)) 99 warning("size is not an integer") 100 if (!missing(size) && length(size) > 1) { # recurse: 101 size = rep(size, length.out = length(x)) 102 ret = lapply(1:length(x), function(i) st_sample(x[i], size[i], type = type, exact = exact, ...)) 103 st_set_crs(do.call(c, ret), st_crs(x)) 104 } else { 105 res = switch(max(st_dimension(x)) + 1, 106 st_multipoints_sample(do.call(c, x), size = size, ..., type = type), 107 st_ll_sample(st_cast(x, "LINESTRING"), size = size, ..., type = type), 108 st_poly_sample(x, size = size, ..., type = type, by_polygon = by_polygon)) 109 if (exact & type == "random" & all(st_geometry_type(res) == "POINT")) { 110 diff = size - length(res) 111 if (diff > 0) { # too few points 112 res_additional = st_sample_exact(x = x, size = diff, ..., 113 type = type, by_polygon = by_polygon) 114 res = c(res, res_additional) 115 } else if (diff < 0) { # too many points 116 res = res[1:size] 117 } 118 } 119 res 120 } 121} 122 123#' @export 124#' @name st_sample 125st_sample.sfg = function(x, size, ...) { 126 st_sample(st_geometry(x), size, ...) 127} 128 129st_poly_sample = function(x, size, ..., type = "random", 130 offset = st_sample(st_as_sfc(st_bbox(x)), 1)[[1]], 131 by_polygon = FALSE) { 132 133 if (by_polygon && inherits(x, "sfc_MULTIPOLYGON")) { # recurse into polygons: 134 sum_a = units::drop_units(sum(st_area(x))) 135 x = lapply(suppressWarnings(st_cast(st_geometry(x), "POLYGON")), st_sfc, crs = st_crs(x)) 136 a = sapply(x, st_area) 137 ret = mapply(st_poly_sample, x, size = size * a / sum_a, type = type, ...) 138 do.call(c, ret) 139 } else if (type %in% c("hexagonal", "regular", "random")) { 140 141 if (isTRUE(st_is_longlat(x))) { 142 if (type == "regular") { 143 message_longlat("st_sample") 144 x = st_set_crs(x, NA) 145 } 146 if (type == "hexagonal") 147 stop("hexagonal sampling on geographic coordinates not supported; consider projecting first") 148 } 149 150 a0 = as.numeric(st_area(st_make_grid(x, n = c(1,1)))) 151 a1 = as.numeric(sum(st_area(x))) 152 # st_polygon(list(rbind(c(-180,-90),c(180,-90),c(180,90),c(-180,90),c(-180,-90)))) 153 # for instance has 0 st_area 154 if (is.finite(a0) && is.finite(a1) && a0 > a0 * 0.0 && a1 > a1 * 0.0) { 155 r = round(size * a0 / a1) 156 size = if (r == 0) 157 rbinom(1, 1, size * a0 / a1) 158 else 159 r 160 } 161 bb = st_bbox(x) 162 163 pts = if (type == "hexagonal") { 164 dx = sqrt(a0 / size / (sqrt(3)/2)) 165 hex_grid_points(x, pt = offset, dx = dx) 166 } else if (type == "regular") { 167 dx = as.numeric(sqrt(a0 / size)) 168 offset = c((offset[1] - bb["xmin"]) %% dx, 169 (offset[2] - bb["ymin"]) %% dx) + bb[c("xmin", "ymin")] 170 n = c(round((bb["xmax"] - offset[1])/dx), round((bb["ymax"] - offset[2])/dx)) 171 st_make_grid(x, cellsize = c(dx, dx), offset = offset, n = n, what = "corners") 172 } else if (type == "random") { 173 lon = runif(size, bb[1], bb[3]) 174 lat = if (isTRUE(st_is_longlat(x))) { # sampling on the sphere: 175 toRad = pi/180 176 lat0 = (sin(bb[2] * toRad) + 1)/2 177 lat1 = (sin(bb[4] * toRad) + 1)/2 178 y = runif(size, lat0, lat1) 179 asin(2 * y - 1) / toRad # http://mathworld.wolfram.com/SpherePointPicking.html 180 } else 181 runif(size, bb[2], bb[4]) 182 m = cbind(lon, lat) 183 st_sfc(lapply(seq_len(nrow(m)), function(i) st_point(m[i,])), crs = st_crs(x)) 184 } 185 pts[x] # cut out x from bbox 186 } else { # try to go into spatstat 187 if (!requireNamespace("spatstat.core", quietly = TRUE)) 188 stop("package spatstat.core required, please install it (or the full spatstat package) first") 189 spatstat_fun = try(get(paste0("r", type), asNamespace("spatstat.core")), silent = TRUE) 190 if (inherits(spatstat_fun, "try-error")) 191 stop(paste0("r", type), " is not an exported function from spatstat.core.") 192 pts = try(spatstat_fun(..., win = spatstat.geom::as.owin(x)), silent = TRUE) 193 if (inherits(pts, "try-error")) 194 stop("The spatstat function ", paste0("r", type), 195 " did not return a valid result. Consult the help file.\n", 196 "Error message from spatstat:\n", pts) 197 st_as_sf(pts)[-1,] 198 } 199} 200 201st_multipoints_sample = function(x, size, ..., type = "random") { 202 if (!inherits(x, "MULTIPOINT")) 203 stop("points sampling only implemented for MULTIPOINT; use sample to sample individual features", call.=FALSE) 204 m = unclass(x) 205 st_sfc(st_multipoint(m[sample(nrow(m), size, ...),]), crs = st_crs(x)) 206} 207 208st_ll_sample = function (x, size, ..., type = "random", offset = runif(1)) { 209 crs = st_crs(x) 210 if (isTRUE(st_is_longlat(x))) { 211 message_longlat("st_sample") 212 st_crs(x) = NA_crs_ 213 } 214 l = st_length(x) 215 if (inherits(l, "units")) 216 l = drop_units(l) 217 if (type == "random") { 218 d = runif(size, 0, sum(l)) 219 } else if (type == "regular") { 220 d = ((1:size) - (1. - (offset %% 1)))/size * sum(l) 221 } else { 222 stop(paste("sampling type", type, "not available for LINESTRING")) # nocov 223 } 224 lcs = c(0, cumsum(l)) 225 if (sum(l) == 0) { 226 grp = list(0) # nocov 227 message("line is of length zero, only one point is sampled") # nocov 228 } else { 229 grp = split(d, cut(d, lcs, include.lowest = TRUE)) 230 grp = lapply(seq_along(x), function(i) grp[[i]] - lcs[i]) 231 } 232 st_sfc(CPL_gdal_linestring_sample(x, grp), crs = crs) 233} 234 235### return points on a triangular grid that 236## - covers a bounding box st_bbox(obj) 237## - contains pt 238## - has x spacing dx: the shortest distance between x coordinates with identical y coordinate 239hex_grid_points = function(obj, pt, dx) { 240 241 bb = st_bbox(obj) 242 dy = sqrt(3) * dx / 2 243 xlim = bb[c("xmin", "xmax")] 244 ylim = bb[c("ymin", "ymax")] 245 offset = c(x = (pt[1] - xlim[1]) %% dx, y = (pt[2] - ylim[1]) %% (2 * dy)) 246 x = seq(xlim[1] - dx, xlim[2] + dx, dx) + offset[1] 247 y = seq(ylim[1] - 2 * dy, ylim[2] + 2 * dy, dy) + offset[2] 248 249 y <- rep(y, each = length(x)) 250 x <- rep(c(x, x + dx / 2), length.out = length(y)) 251 xy = cbind(x, y)[x >= xlim[1] & x <= xlim[2] & y >= ylim[1] & y <= ylim[2], ] 252 st_sfc(lapply(seq_len(nrow(xy)), function(i) st_point(xy[i,])), crs = st_crs(bb)) 253} 254 255st_sample_exact = function(x, size, ..., type, by_polygon) { 256 random_pt = st_sample(x = x, size = size, ..., type = type, exact = FALSE) 257 while (length(random_pt) < size) { 258 diff = size - length(random_pt) 259 random_pt_new = st_sample(x, size = diff, ..., type, exact = FALSE, by_polygon = by_polygon) 260 random_pt = c(random_pt, random_pt_new) 261 } 262 if(length(random_pt) > size) { 263 random_pt = random_pt[1:size] 264 } 265 random_pt 266} 267