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