1# Author: Robert J. Hijmans
2# Date : November 2018
3# Version 1.0
4# License GPL v3
5
6
7.big_number_warning <- function() {
8# this warning should be given by C
9	warn("big number", "cell numbers larger than ", 2^.Machine$double.digits, " are approximate")
10}
11
12
13ext_from_rc <- function(x, r1, r2, c1, c2){
14	e <- as.vector(ext(x))
15	r <- res(x)
16	c1 <- min(max(c1, 1), ncol(x))
17	c2 <- min(max(c2, 1), ncol(x))
18	if (c1 > c2) {
19		tmp <- c1
20		c1 <- c2
21		c2 <- tmp
22	}
23	r1 <- min(max(r1, 1), nrow(x))
24	r2 <- min(max(r2, 1), nrow(x))
25	if (r1 > r2) {
26		tmp <- r1
27		r1 <- r2
28		r2 <- tmp
29	}
30
31	xn <- xFromCol(x, c1) - 0.5 * r[1]
32	xx <- xFromCol(x, c2) + 0.5 * r[1]
33	yx <- yFromRow(x, r1) + 0.5 * r[2]
34	yn <- yFromRow(x, r2) - 0.5 * r[2]
35	ext(c(sort(c(xn, xx))), sort(c(yn, yx)))
36}
37
38
39
40.makeDataFrame <- function(x, v, factors=TRUE) {
41	v <- data.frame(v, check.names = FALSE)
42	if (factors) {
43		ff <- is.factor(x)
44		if (any(ff)) {
45			ff <- which(ff)
46			#levs <- levels(x)
47			cgs <- cats(x)
48			for (f in ff) {
49				#lvs <- levs[[f]]
50				cg <- cgs[[f]]
51				i <- match(v[,f], cg[,1])
52				act <- activeCat(x, f) + 1
53				#v[[f]] = factor(v[[f]], levels=(1:length(lvs))-1)
54				#levels(v[[f]]) = levs[[f]]
55
56				if (!inherits(cg[[act]], "numeric")) {
57					v[[f]] <- factor(cg[i, act], levels=unique(cg[[act]]))
58				} else {
59					v[[f]] <- cg[i, act]
60				}
61			}
62		}
63	}
64	v
65}
66
67setlabs <- function(x, labs) {
68	x[ (x<1) | (x>length(labs))] <- NA
69	x <- factor(x, levels=1:length(labs))
70	levels(x) <- labs
71	x
72}
73
74
75wmean <- function(p, na.rm=FALSE) {
76	n <- length(p)
77	w <- p[[n]]
78	p[[n]] <- NULL
79	sapply(p, function(x) {
80		stats::weighted.mean(x, w, na.rm=na.rm)
81	})
82}
83
84wsum <- function(p, na.rm=FALSE) {
85	n <- length(p)
86	w <- p[[n]]
87	p[[n]] <- NULL
88	sapply(p, function(x) {
89		sum(x * w, na.rm=na.rm)
90	})
91}
92
93wmin <- function(p, na.rm=FALSE) {
94	n <- length(p)
95	p[[n]] <- NULL
96	sapply(p, function(x) {
97		min(x, na.rm=na.rm)
98	})
99}
100
101wmax <- function(p, na.rm=FALSE) {
102	n <- length(p)
103	p[[n]] <- NULL
104	sapply(p, function(x) {
105		max(x, na.rm=na.rm)
106	})
107}
108
109
110		#if (!list) {
111			#if (geomtype(y) == "points")  {
112			#	e <- cbind(ID=1:length(e), matrix(unlist(e), ncol=nlyr(x), byrow=TRUE))
113			#} else {
114			#	e <- lapply(1:length(e), function(i) {
115			#		ee <- unlist(e[[i]])
116			#		if (length(ee) == 0) ee <- NA
117			#		cbind(ID=i, matrix(ee, ncol=length(e[[i]])))
118			#	})
119			#	e <- do.call(rbind, e)
120			#}
121		#}
122
123
124
125setMethod("extract", signature(x="SpatRaster", y="SpatVector"),
126function(x, y, fun=NULL, method="simple", list=FALSE, factors=TRUE, cells=FALSE, xy=FALSE, weights=FALSE, exact=FALSE, touches=is.lines(y), layer=NULL, ...) {
127
128	nl <- nlyr(x)
129	useLyr <- FALSE
130	method <- match.arg(tolower(method), c("simple", "bilinear"))
131	hasfun <- !is.null(fun)
132	if (weights && exact) {
133		exact = FALSE
134	}
135	if (hasfun) {
136		cells <- FALSE
137		xy <- FALSE
138		if (weights || exact) {
139			list <- TRUE
140			fun <- .makeTextFun(fun)
141			bad <- FALSE
142			if (is.character(fun)) {
143				if (!(fun %in% c("sum", "mean", "min", "max"))) {
144					bad <- TRUE
145				} else if (fun == "mean") {
146					fun <- wmean
147				} else if (fun == "sum") {
148					fun <- wsum
149				} else if (fun == "min") {
150					fun <- wmin
151				} else if (fun == "max") {
152					fun <- wmax
153				}
154			} else {
155				bad <- TRUE
156			}
157			if (bad) error("extract", 'if weights or exact=TRUE, "fun" must be "sum", "mean" or NULL')
158		}
159	}
160	if (!is.null(layer) && nl > 1) {
161		if (any(is.na(layer))) {error("extract", "argument 'layer' cannot have NAs")}
162		if (length(layer) == 1) {
163			lyr_name <- layer
164			layer <- as.character(y[[layer,drop=TRUE]])
165		} else {
166			lyr_name <- "layer"
167			stopifnot(length(layer) == nrow(y))
168		}
169		if (is.numeric(layer)) {
170			layer <- round(layer)
171			stopifnot(min(layer) > 0 & max(layer) <= nl)
172		} else {
173			layer <- match(layer, names(x))
174			if (any(is.na(layer))) error("extract", "names in argument 'layer' do not match names(x)")
175		}
176		useLyr <- TRUE
177	}
178
179	#f <- function(i) if(length(i)==0) { NA } else { i }
180	#e <- rapply(e, f, how="replace")
181	cn <- names(x)
182	opt <- spatOptions()
183	if (list) {
184		e <- x@ptr$extractVector(y@ptr, touches[1], method, isTRUE(cells[1]), isTRUE(xy[1]), isTRUE(weights[1]), isTRUE(exact[1]), opt)
185		x <- messages(x, "extract")
186		if (weights || exact) {
187			if (hasfun) {
188				e <- sapply(e, fun, ...)
189				e <- matrix(e, nrow=nrow(y), byrow=TRUE)
190				colnames(e) <- cn
191				e <- cbind(ID=1:nrow(e), e)
192			}
193		}
194		return(e)
195	}
196
197	e <- x@ptr$extractVectorFlat(y@ptr, touches[1], method, isTRUE(cells[1]), isTRUE(xy[1]), isTRUE(weights[1]), isTRUE(exact[1]), opt)
198	x <- messages(x, "extract")
199	nc <- nl
200	if (cells) {
201		cn <- c(cn, "cell")
202		nc <- nc + 1
203	}
204	if (weights) {
205		cn <- c(cn, "weight")
206		nc <- nc + 1
207	} else if (exact) {
208		cn <- c(cn, "fraction")
209		nc <- nc + 1
210	}
211	if (xy) {
212		cn <- c(cn, "x", "y")
213		nc <- nc + 2
214	}
215
216	geo <- geomtype(y)
217	if (geo == "points") {
218		## this was? should be fixed upstream
219		if (nc == nl) {
220			e <- matrix(e, ncol=nc)
221		} else {
222			e <- matrix(e, ncol=nc, byrow=TRUE)
223		}
224		e <- cbind(1:nrow(e), e)
225		if (nrow(e) > nrow(y)) { #multipoint
226			g <- geom(y)
227			e[,1] <- g[,1]
228		}
229	} else {
230		e <- matrix(e, ncol=nc+1, byrow=TRUE)
231	}
232	cn <- c("ID", cn)
233	colnames(e) <- cn
234	if (hasfun) {
235		fun <- match.fun(fun)
236		e <- data.frame(e)
237		e <- aggregate(e[,-1,drop=FALSE], e[,1,drop=FALSE], fun, ...)
238
239		m <- sapply(e, NCOL)
240		if (any(m > 1)) {
241			e <- do.call(cbind, as.list(e))
242			skip <- (length(cn) - nlyr(x))
243			nms <- colnames(e)
244			snms <- nms[(skip+1):length(nms)]
245			mr <- max(m)
246			if (!all(snms=="")) {
247				snms <- paste0(rep(names(x), each=mr), ".", snms)
248			} else {
249				snms <- paste0(rep(names(x), each=mr), ".", rep(1:mr))
250			}
251			snms <- c(cn[1:skip], snms)
252			colnames(e) <- snms
253			e <- data.frame(e)
254		}
255	} else if (cells) {
256		cncell <- cn =="cell"
257		e[, cncell] <- e[, cncell] + 1
258	}
259
260	if (factors) {
261		id <- data.frame(e[,1,drop=FALSE])
262		e <- cbind(id, .makeDataFrame(x, e[,-1,drop=FALSE], TRUE))
263	}
264
265	if (useLyr) {
266		idx <- cbind(e[,1], layer[e[,1]]+1)
267		ee <- cbind(e[,1,drop=FALSE], names(x)[idx[,2]-1], value=e[idx])
268		colnames(ee)[2] <- lyr_name
269		if (ncol(e) > (nl+1)) {
270			e <- cbind(ee, e[,(nl+1):ncol(e), drop=FALSE])
271		} else {
272			e <- ee
273		}
274	}
275	e
276})
277
278
279setMethod("[", c("SpatRaster", "SpatVector", "missing"),
280function(x, i, j, ... , drop=FALSE) {
281	v <- extract(x, i)
282	if (drop) {
283		as.vector(v)
284	} else {
285		v
286	}
287})
288
289setMethod("[", c("SpatVector", "SpatVector", "missing"),
290function(x, i, j, ... , drop=FALSE) {
291	r <- !relate(x, i, "disjoint")
292	r <- which(apply(r, 1, any))
293	x[r, ]
294})
295
296
297setMethod("[", c("SpatVector", "SpatExtent", "missing"),
298function(x, i, j, ... , drop=FALSE) {
299	x[as.polygons(i)]
300})
301
302
303setMethod("extract", signature(x="SpatRaster", y="matrix"),
304function(x, y, ...) {
305	.checkXYnames(colnames(y))
306	#if (length(list(...)) == 0) {
307	#	i <- cellFromXY(x, y)
308	#	r <- cbind(1:length(i), x[i])
309	#	colnames(r) <- c("ID", names(x))
310	#} else {
311	y <- vect(y)
312	e <- extract(x, y, ...)
313	if (NCOL(e) > 1) {
314		e[, -1, drop=FALSE]
315	} else {
316		e
317	}
318	#}
319})
320
321
322setMethod("extract", signature(x="SpatRaster", y="data.frame"),
323function(x, y, ...) {
324	if (ncol(y) != 2) {
325		error("extract", "extract expects a 2 column data.frame of x and y coordinates")
326	}
327	v <- vect(y, colnames(y))
328	extract(x, v, ...)
329})
330
331
332setMethod("extract", signature(x="SpatRaster", y="numeric"),
333function(x, y, ...) {
334	y <- as.integer(y)
335	y[(y < 1) | (y > ncell(x))] <- NA
336	x[y]
337})
338
339setMethod("extract", signature(x="SpatRaster", y="SpatExtent"),
340function(x, y, factors=TRUE, cells=FALSE, xy=FALSE) {
341	y <- cells(x, y)
342	if (factors) dataframe = TRUE
343	v <- extract_cell(x, y, factors=factors)
344	if (cells) {
345		v$cell <- y
346	}
347	if (xy) {
348		v <- cbind(v, xyFromCell(x, y))
349	}
350	v
351}
352)
353
354
355setMethod("[", c("SpatRaster", "missing", "missing"),
356function(x, i, j, ... , drop=FALSE) {
357	values(x, mat=!drop)
358})
359
360setMethod("[", c("SpatRaster", "logical", "missing"),
361function(x, i, j, ... , drop=FALSE) {
362	x[which(i),, drop=drop]
363})
364
365
366extract_cell <- function(x, cells, drop=FALSE, factors=TRUE) {
367	e <- x@ptr$extractCell(cells-1)
368	messages(x, "extract_cell")
369	e <- do.call(cbind, e)
370	colnames(e) <- names(x)
371	.makeDataFrame(x, e, factors)[,,drop]
372}
373
374
375setMethod("[", c("SpatRaster", "numeric", "missing"),
376function(x, i, j, ... , drop=TRUE) {
377
378	add <- any(grepl("drop", names(match.call())))
379	if (!drop) {
380		if (nargs() == 3) {
381			rc <- rowColFromCell(x, i)
382			e <- ext_from_rc(x, min(rc[,1]), max(rc[,1]), min(rc[,2]), max(rc[,2]))
383		} else {
384			e <- ext_from_rc(x, min(i), max(i), 1, ncol(x))
385		}
386		return(crop(x, e))
387	}
388	if (nargs() > (2+add)) {
389		i <- cellFromRowColCombine(x, i, 1:ncol(x))
390	}
391	extract_cell(x, i, drop=FALSE)
392})
393
394
395setMethod("[", c("SpatRaster", "missing", "numeric"),
396function(x, i, j, ... , drop=TRUE) {
397	if (!drop) {
398		e <- ext_from_rc(x, 1, nrow(x), min(j), max(j))
399		return(crop(x, e))
400	}
401
402	i <- cellFromRowColCombine(x, 1:nrow(x), j)
403	extract_cell(x, i, drop=FALSE)
404})
405
406
407setMethod("[", c("SpatRaster", "numeric", "numeric"),
408function(x, i, j, ..., drop=TRUE) {
409	if (!drop) {
410		e <- ext_from_rc(x, min(i), max(i), min(j), max(j))
411		return(crop(x, e))
412	}
413	i <- cellFromRowColCombine(x, i, j)
414	extract_cell(x, i, drop=FALSE)
415})
416
417
418setMethod("[", c("SpatRaster", "SpatRaster", "missing"),
419function(x, i, j, ..., drop=TRUE) {
420	x <- x[ext(i), drop=drop]
421	if (!drop) {
422		if (compareGeom(x, i, crs=FALSE, stopOnError=FALSE)) {
423			x <- mask(x, i)
424		}
425	}
426	x
427})
428
429setMethod("[", c("SpatRaster", "SpatExtent", "missing"),
430function(x, i, j, ..., drop=FALSE) {
431	x <- crop(x, i)
432	if (drop) {
433		values(x)
434	} else {
435		x
436	}
437})
438
439
440
441setMethod("extract", c("SpatVector", "SpatVector"),
442function(x, y, ...) {
443	r <- relate(y, x, "within")
444	e <- apply(r, 1, which)
445	if (length(e) == 0) {
446		e <- list(e)
447	}
448	if (is.list(e)) {
449		e <- lapply(1:length(e), function(i) {
450			if (length(e[[i]]) == 0) {
451				cbind(i, NA)
452			} else {
453				cbind(i, e[[i]])
454			}
455		})
456		e <- do.call(rbind, e)
457	} else {
458		e <- cbind(1:nrow(y), e)
459	}
460	if (ncol(x) > 0) {
461		d <- as.data.frame(x)
462		e <- data.frame(id.x=e[,1], d[e[,2], ,drop=FALSE])
463		rownames(e) <- NULL
464	} else {
465		colnames(e) <- c("id.x", "id.y")
466	}
467	e
468})
469
470