1 2.seed <- function() { 3 sample.int(.Machine$integer.max, 1) 4} 5 6 7.sampleCells <- function(x, size, method, replace, na.rm=FALSE, ext=NULL) { 8 r <- rast(x) 9 lonlat <- is.lonlat(r, perhaps=TRUE, warn=TRUE) 10 if (!is.null(ext)) { 11 r <- crop(rast(r), ext) 12 } 13 if ( ((!replace) || (method == "regular")) && (size >= ncell(r)) ) { 14 cells <- 1:ncell(r) 15 } 16 17 if (method == "random") { 18 nsize <- size 19 if (na.rm) { 20 if (replace) { 21 size <- size*5 22 } else { 23 size <- min(ncell(r)*2, size*5) 24 } 25 } 26 if (lonlat) { 27 m <- ifelse(replace, 1.5, 1.25) 28 n <- m * size 29 y <- yFromRow(r, 1:nrow(r)) 30 w <- abs(cos(pi*y/180)) 31 rows <- sample.int(nrow(r), n, replace=TRUE, prob=w) 32 cols <- sample.int(ncol(r), n, replace=TRUE) 33 cells <- cellFromRowCol(r, rows, cols) 34 if (!replace) { 35 cells <- unique(cells) 36 } 37 } else { 38 cells <- sample(ncell(r), size, replace=replace) 39 } 40 } else { # regular 41 if (lonlat) { 42 ratio <- 0.5 * ncol(r)/nrow(r) 43 n <- sqrt(size) 44 nx <- max(1, (round(n*ratio))) 45 ny <- max(1, (round(n/ratio))) 46 xi <- ncol(r) / nx 47 yi <- nrow(r) / ny 48 rows <- unique(round(seq(.5*yi, nrow(r), yi))) 49 50 w <- cos(pi*yFromRow(r, rows)/180) 51 w <- w * length(w)/sum(w) 52 xi <- xi / w 53 xi <- pmax(1,pmin(xi, ncol(r))) 54 z <- list() 55 #off <- stats::runif(1) 56 for (i in 1:length(rows)) { 57 z[[i]] <- cbind(rows[i], unique(round(seq(0.5*xi[i], ncol(r), xi[i])))) 58 } 59 z <- do.call(rbind, z) 60 cells <- cellFromRowCol(r, z[,1], z[,2]) 61 62 } else { 63 f <- sqrt(size / ncell(r)) 64 nr <- ceiling(nrow(r) * f) 65 nc <- ceiling(ncol(r) * f); 66 xstep <- ncol(r) / nc 67 ystep <- nrow(r) / nr 68 xsamp <- seq(0.5*xstep, ncol(r), xstep) 69 ysamp <- seq(0.5*ystep, nrow(r), ystep) 70 xy <- expand.grid(ysamp, xsamp) 71 cells <- cellFromRowCol(r, xy[,1], xy[,2]) 72 } 73 } 74 if (!is.null(ext)) { 75 cells <- cellFromXY(x, xyFromCell(r, cells)) 76 } 77 if (na.rm) { 78 v <- rowSums(is.na(x[cells])) == 0 79 cells <- cells[v] 80 } 81 if (method == "random") { 82 if (length(cells) > nsize) { 83 cells <- cells[1:nsize] 84 } 85 } 86 return(cells) 87} 88 89 90set_factors <- function(x, ff, cts, asdf) { 91 if (any(ff)) { 92 x <- data.frame(x) 93 for (i in which(ff)) { 94 ct <- cts[[i]] 95 m <- match(x[[i]], ct[,1]) 96 if (!inherits(ct[[2]], "numeric")) { 97 x[[i]] <- factor(ct[m,2], levels=unique(ct[[2]])) 98 } else { 99 x[[i]] <- ct[m,2] 100 } 101 } 102 } else if (asdf) { 103 x <- data.frame(x) 104 } 105 x 106} 107 108 109setMethod("spatSample", signature(x="SpatRaster"), 110 function(x, size, method="random", replace=FALSE, na.rm=FALSE, as.raster=FALSE, as.df=TRUE, as.points=FALSE, values=TRUE, cells=FALSE, xy=FALSE, ext=NULL, warn=TRUE) { 111 112 size <- round(size) 113 if (any(size < 1)) { 114 error("spatSample", "sample size must be a positive integer") 115 } 116 if ((size > ncell(x)) & (!replace)) { 117 size <- ncell(x) 118 } 119 120 if (!as.raster) { 121 ff <- is.factor(x) 122 lv <- active_cats(x) 123 } 124 125 if (cells || xy || as.points) { 126 size <- size[1] 127 cnrs <- .sampleCells(x, size, method, replace, na.rm, ext) 128 if (method == "random") { 129 if (length(cnrs) < size) { 130 warn("spatSample", "fewer cells returned than requested") 131 } else if (length(cnrs) > size) { 132 cnrs <- cnrs[1:size] 133 } 134 } 135 out <- NULL 136 if (cells) { 137 out <- matrix(cnrs, ncol=1) 138 colnames(out) <- "cell" 139 } 140 if (xy) { 141 out <- cbind(out, xyFromCell(x, cnrs)) 142 } 143 if (values && hasValues(x)) { 144 e <- extract(x, cnrs) 145 e <- set_factors(e, ff, lv, as.df) 146 if (is.null(out)) { 147 out <- e 148 } else { 149 out <- cbind(out, e) 150 } 151 } 152 if (as.points) { 153 if (xy) { 154 out <- vect(out, crs=crs(x)) 155 } else { 156 xy <- xyFromCell(x, cnrs) 157 # xy is a matrix, no geom argument 158 v <- vect(xy, crs=crs(x)) 159 values(v) <- out 160 return(v) 161 } 162 } 163 return(out) 164 } 165 if (!hasValues(x) & !as.raster) { 166 error("spatSample", "SpatRaster has no values") 167 } 168 169 method <- tolower(method) 170 stopifnot(method %in% c("random", "regular")) 171 if (!replace) size <- pmin(ncell(x), size) 172 173 if (!is.null(ext)) x <- crop(x, ext) 174 175 if (method == "regular") { 176 if (as.raster) { 177 if (length(size) > 1) { 178 x@ptr <- x@ptr$sampleRowColRaster(size[1], size[2]) 179 } else { 180 x@ptr <- x@ptr$sampleRegularRaster(size) 181 } 182 x <- messages(x, "spatSample") 183 return(x); 184 } else { 185 opt <- spatOptions() 186 if (length(size) > 1) { 187 v <- x@ptr$sampleRowColValues(size[1], size[2], opt) 188 } else { 189 v <- x@ptr$sampleRegularValues(size, opt) 190 } 191 x <- messages(x, "spatSample") 192 if (length(v) > 0) { 193 v <- do.call(cbind, v) 194 colnames(v) <- names(x) 195 } 196 v <- set_factors(v, ff, lv, as.df) 197 return(v) 198 } 199 } else { # random 200 size <- size[1] 201 if (as.raster) { 202 x@ptr <- x@ptr$sampleRandomRaster(size, replace, .seed()) 203 x <- messages(x, "spatSample") 204 return(x); 205 } else { 206 #v <- x@ptr$sampleRandomValues(size, replace, seed) 207 if (size > 0.75 * ncell(x)) { 208 if (na.rm) { 209 out <- stats::na.omit(values(x)) 210 attr(x, "na.action") <- NULL 211 if (nrow(out) < size) { 212 warn("spatSample", "more non-NA cells requested than available") 213 } else { 214 out <- out[sample(nrow(out), size), ,drop=FALSE] 215 } 216 } else { 217 out <- values(x) 218 out <- out[sample(nrow(out), size, replace=replace), ,drop=FALSE] 219 } 220 out <- set_factors(out, ff, lv, as.df) 221 return(out) 222 } 223 224 if (na.rm) { 225 scells <- NULL 226 ssize <- size*2 227 for (i in 1:10) { 228 scells <- c(scells, .sampleCells(x, ssize, method, replace)) 229 if ((i>1) && (!replace)) { 230 scells <- unique(scells) 231 } 232 out <- stats::na.omit(x[scells]) 233 if (nrow(out) >= size) { 234 out <- out[1:size, ,drop=FALSE] 235 attr(out, "na.action") <- NULL 236 rownames(out) <- NULL 237 break 238 } 239 } 240 } else { 241 scells <- .sampleCells(x, size, method, replace) 242 out <- x[scells] 243 } 244 if (NROW(out) < size) { 245 if (warn) warn("spatSample", "fewer values returned than requested") 246 } else if (method == "random") { 247 if (is.null(dim(out))) { 248 out = out[1:size] 249 } else { 250 out = out[1:size, ,drop=FALSE] 251 } 252 } 253 return(out) 254 } 255 } 256 } 257) 258 259 260setMethod("spatSample", signature(x="SpatExtent"), 261 function(x, size, method="random", lonlat, as.points=FALSE) { 262 if (missing(lonlat)) { 263 error("spatSample", "provide a lonlat argument") 264 } 265 if (lonlat) { 266 stopifnot(x$ymax <= 90 || x$ymin >= -90) 267 } 268 method <- match.arg(method, c("regular", "random")) 269 size <- round(size) 270 stopifnot(size > 0) 271 if (method=="random") { 272 s <- x@ptr$sampleRandom(size, lonlat, .seed()) 273 } else { 274 s <- x@ptr$sampleRegular(size, lonlat) 275 } 276 s <- do.call(cbind, s) 277 colnames(s) <- c("x", "y") 278 if (as.points) { 279 s <- vect(s) 280 } 281 s 282 } 283) 284 285 286 287 288 289.grid_sample <- function(xy, n=1, r, chess="") { 290 291 cell <- cellFromXY(r, xy) 292 uc <- unique(stats::na.omit(cell)) 293 294 chess <- trim(chess) 295 if (chess != "") { 296 chess <- match.arg(tolower(chess), c("white", "black")) 297 nc <- ncol(r) 298 if (nc %% 2 == 1) { 299 if (chess=="white") { 300 tf <- 1:ceiling(ncell(r)/2) * 2 - 1 301 } else { 302 tf <- 1:ceiling((ncell(r)-1)/2) * 2 303 } 304 } else { 305 nr <- nrow(r) 306 row1 <- 1:(ceiling(nr / 2)) * 2 - 1 307 row2 <- row1 + 1 308 row2 <- row2[row2 <= nr] 309 310 if (chess=="white") { 311 col1 <- 1:(ceiling(nc / 2)) * 2 - 1 312 col2 <- col1 + 1 313 col2 <- col2[col2 <= nc] 314 } else { 315 col1 <- 1:(ceiling(nc / 2)) * 2 316 col2 <- col1 - 1 317 col1 <- col1[col1 <= nc] 318 } 319 320 cells1 <- cellFromRowColCombine(r, row1, col1) 321 cells2 <- cellFromRowColCombine(r, row2, col2) 322 tf <- c(cells1, cells2) 323 } 324 uc <- uc[uc %in% tf] 325 } 326 327 cell <- cellFromXY(r, xy) 328 cell <- cbind(1:nrow(xy), cell, stats::runif(nrow(xy))) 329 cell <- stats::na.omit(cell) 330 331 cell <- cell[order(cell[,3]), ] 332 sel <- list() 333 for (i in 1:length(uc)) { 334 ss <- subset(cell, cell[,2] == uc[i]) 335 sel[[i]] <- ss[1:min(n, nrow(ss)), 1] 336 } 337 unlist(sel) 338} 339 340 341#coordinates <- function(x) { 342# do.call(cbind, x@ptr$coordinates()) 343#} 344 345get_field_name <- function(x, nms, sender="") { 346 x <- x[1] 347 if (is.numeric(x)) { 348 x <- round(x) 349 if (x > 0 && x <= length(nms)) { 350 x = nms[x] 351 } else { 352 error(sender, "invalid index. there are ", length(nms), " columns") 353 } 354 } else if (is.character(x)) { 355 if (!(x %in% nms)) { 356 error(sender, "invalid name") 357 } 358 } 359 x 360} 361 362 363setMethod("spatSample", signature(x="SpatVector"), 364 function(x, size, method="random", strata=NULL, chess="") { 365 method = match.arg(tolower(method), c("regular", "random")) 366 stopifnot(size > 0) 367 gtype <- geomtype(x) 368 if (gtype == "polygons") { 369 if (!is.null(strata)) { 370 if (length(strata) == 1) { 371 if (is.character(strata)) { 372 stopifnot(strata %in% names(x)) 373 } else { 374 stopifnot((strata > 0) && (strata < ncol(x))) 375 } 376 strata <- x[[strata, drop=TRUE]] 377 } else if (length(strata) != length(x)) { 378 stop("length of strata must be 1 or length(x)") 379 } 380 s <- stats::na.omit(unique(strata)) 381 n <- length(size) 382 if (n==1) { 383 n <- rep_len(n, length(s)) 384 } else if (length(s) != n) { 385 stop("length of strata must be 1 or length(na.omit(unique(strata)))") 386 } 387 r <- lapply(s, function(s) { 388 spatSample(x[strata == s, ], size, method, NULL, "") 389 }) 390 r <- do.call(rbind, r) 391 return(r) 392 } 393 if (length(size) == 1) { 394 x@ptr = x@ptr$sample(size, method[1], .seed()) 395 } else { 396 x@ptr = x@ptr$sampleGeom(size, method[1], .seed()) 397 } 398 return(messages(x)) 399 } else if (grepl(gtype, "points")) { 400 if (!is.null(strata)) { 401 if (inherits(strata, "SpatRaster")) { 402 xy <- crds(x) 403 i <- .grid_sample(xy, size[1], rast(strata), chess) 404 return(x[i,]) 405 } else { 406 error("spatSample", "not yet implemented for these strata") 407 } 408 } else { 409 error("spatSample", "use `sample` to sample (point) geometries") 410 } 411 } else { 412 error("spatSample", "not yet implemented for lines") 413 } 414 } 415) 416 417#spatSample(disagg(as.points(v)), 1, "stratified", strata=r, chess="") 418 419 420 421# setMethod("spatSample", signature(x="SpatExtent"), 422 # function(x, size, method="regular", lonlat, ...) { 423 # if (missing(lonlat)) { 424 # stop("provide a lonlat argument") 425 # } 426 # method = match.arg(method, c("regular", "random")) 427 # size <- round(size) 428 # stopifnot(size > 0) 429 # e <- as.vector(x) 430 # if (method=="random") { 431 # if (lonlat) { 432 # d <- round((e[4] - e[3]) * 1000); 433 # dx <- (e[4] - e[3]) / (2 * d) 434 # r <- unique(seq(e[3], e[4], length.out=d)) 435 # w <- abs(cos(pi*r/180)) 436 # x <- sample.int(length(r), size, prob=w, replace=TRUE) 437 # lat <- r[x] + stats::runif(size, -dx, dx) 438 # lon <- stats::runif(size, min = e[1], max = e[2]) 439 # vect(cbind(lon,lat), crs="+proj=lonlat +datum=WGS84") 440 # } else { 441 # x <- stats::runif(size, min = e[1], max = e[2]) 442 # y <- stats::runif(size, min = e[3], max = e[4]) 443 # vect(cbind(x, y)) 444 # } 445 # } else { 446 # r <- range(x) 447 # ratio <- 0.5 * r[1]/r[2] 448 # n <- sqrt(size) 449 # nx <- max(1, (round(n*ratio))) 450 # ny <- max(1, (round(n/ratio))) 451 # xi <- r[1] / nx 452 # yi <- r[2] / ny 453 # if (lonlat) { 454 # lat <- seq(e[3]+0.5*yi, e[4], yi) 455 # w <- cos(pi*lat/180) 456 # w <- w * length(w)/sum(w) 457 # xi <- xi / w 458 # xi <- pmin(xi, 180) 459 # z <- list() 460 # #off <- stats::runif(1) 461 # for (i in 1:length(lat)) { 462 # z[[i]] <- cbind(seq(e[1]+0.5*xi[i], e[2], xi[i]), lat[i]) 463 # } 464 # z <- do.call(rbind, z) 465 # vect(z, crs="+proj=lonlat +datum=WGS84") 466 # } else { 467 # x <- seq(e[1]+0.5*xi, e[2], xi) 468 # y <- seq(e[3]+0.5*yi, e[4], yi) 469 # vect(cbind(rep(x, length(y)), rep(y, each=length(x)))) 470 # } 471 # } 472 # } 473# ) 474