1# Author: Robert J. Hijmans
2# Date : December 2009
3# Version 0.9
4# Licence GPL v3
5
6
7
8setMethod('extract', signature(x='Raster', y='SpatialPolygons'),
9function(x, y, fun=NULL, na.rm=FALSE, exact=FALSE, weights=FALSE, normalizeWeights=TRUE, cellnumbers=FALSE, small=TRUE, df=FALSE, layer, nl, factors=FALSE, sp=FALSE, ...){
10
11	#px <-.getCRS(x, asText=FALSE)
12	px <-.getCRS(x)
13	comp <- compareCRS(px,.getCRS(y), unknown=TRUE)
14	if (!comp) {
15		.requireRgdal()
16		warning('Transforming SpatialPolygons to the crs of the Raster')
17		y <- sp::spTransform(y, px)
18	}
19
20	spbb <- sp::bbox(y)
21	rsbb <- bbox(x)
22	addres <- max(res(x))
23	npol <- length(y@polygons)
24	res <- list()
25	res[[npol+1]] <- NA
26
27	if (!is.null(fun)) {
28		cellnumbers <- FALSE
29	    if (weights || exact) {
30			if (!is.null(fun)) {
31				fun <- match.fun(fun)
32				test <- try(methods::slot(fun, 'generic') == 'mean', silent=TRUE)
33				if (!isTRUE(test)) {
34					warning('"fun" was changed to "mean"; other functions cannot be used when "weights=TRUE"' )
35				}
36			}
37			fun <- function(x, ...) {
38				# some complexity here because different layers could
39				# have different NA cells
40				if ( is.null(x) ) {
41					return(rep(NA, nl))
42				}
43				w <- x[,nl+1]
44				x <- x[,-(nl+1), drop=FALSE]
45				x <- x * w
46				w <- matrix(rep(w, nl), ncol=nl)
47				w[is.na(x)] <- NA
48				w <- colSums(w, na.rm=TRUE)
49				x <- apply(x, 1, function(X) { X / w } )
50				if (!is.null(dim(x))) {
51					rowSums(x, na.rm=na.rm)
52				} else {
53					sum(x, na.rm=na.rm)
54				}
55			}
56		}
57
58		if (sp) {
59			df <- TRUE
60		}
61
62		doFun <- TRUE
63
64	} else {
65		if (sp) {
66			sp <- FALSE
67			df <- FALSE
68			warning('argument sp=TRUE is ignored if fun=NULL')
69		#} else if (df) {
70		#	df <- FALSE
71		#	warning('argument df=TRUE is ignored if fun=NULL')
72		}
73
74		doFun <- FALSE
75	}
76
77	if (missing(layer)) {
78		layer <- 1
79	} else {
80		layer <- max(min(nlayers(x), layer), 1)
81	}
82	if (missing(nl)) {
83		nl <- nlayers(x) - layer + 1
84	} else {
85		nl <- max(min(nlayers(x)-layer+1, nl), 1)
86	}
87
88
89	if (spbb[1,1] >= rsbb[1,2] | spbb[1,2] <= rsbb[1,1] | spbb[2,1] >= rsbb[2,2] | spbb[2,2] <= rsbb[2,1]) {
90		if (df) {
91			res <- data.frame(matrix(ncol=1, nrow=0))
92			colnames(res) <- 'ID'
93			return(res)
94		}
95		return(res[1:npol])
96	}
97
98
99
100	rr <- raster(x)
101
102	pb <- pbCreate(npol, label='extract', ...)
103
104	if (.doCluster()) {
105		cl <- getCluster()
106		on.exit( returnCluster() )
107		nodes <- min(npol, length(cl))
108		message('Using cluster with ', nodes, ' nodes')
109		utils::flush.console()
110
111		.sendCall <- eval( parse( text="parallel:::sendCall") )
112		parallel::clusterExport(cl, c('rsbb', 'rr', 'weights', 'exact', 'addres', 'cellnumbers', 'small'), envir=environment())
113
114		clFun <- function(i, pp) {
115			spbb <- sp::bbox(pp)
116
117			if (spbb[1,1] >= rsbb[1,2] | spbb[1,2] <= rsbb[1,1] | spbb[2,1] >= rsbb[2,2] | spbb[2,2] <= rsbb[2,1]) {
118				# do nothing; res[[i]] <- NULL
119			} else {
120				rc <- crop(rr, extent(pp)+addres)
121				if (weights) {
122					rc <- .polygonsToRaster(pp, rc, getCover=TRUE, silent=TRUE)
123					rc[rc==0] <- NA
124					xy <- rasterToPoints(rc)
125					if (normalizeWeights) {
126						weight <- xy[,3] / sum(xy[,3])
127					} else {
128						weight <- xy[,3] #/ 100
129					}
130					xy <- xy[, -3, drop=FALSE]
131				} else if (exact) {
132					erc <- crop(x, rc)
133					xy <- exactextractr::exact_extract(erc, pp, include_cell=cellnumbers, progress=FALSE)[[1]]
134				} else {
135					rc <- .polygonsToRaster(pp, rc, silent=TRUE)
136					r <- rasterToPoints(rc)[,-3,drop=FALSE]
137				}
138
139				if (length(xy) > 0) { # catch very small polygons
140					if (exact) {
141						if (weights) {
142							if (normalizeWeights) {
143								xy$coverage_fraction <- xy$coverage_fraction / sum(xy$coverage_fraction)
144							}
145							colnames(xy)[ncol(xy)] <- "weight"
146						} else {
147							xy$coverage_fraction  <- NULL
148						}
149						if (cellnumbers) {
150							nms <- colnames(xy)
151							# not good if there is a layer called cell
152							nms <- c("cell", nms[nms != "cell"])
153							xy <- xy[,nms]
154						}
155						r <- as.matrix(xy)
156
157					} else {
158						r <- .xyValues(x, xy, layer=layer, nl=nl)
159						if (weights) {
160							if (cellnumbers) {
161								cell <- cellFromXY(x, xy)
162								r <- cbind(cell, r, weight)
163							} else {
164								r <- cbind(r, weight)
165							}
166						} else if (cellnumbers) {
167							cell <- cellFromXY(x, xy)
168							r <- cbind(cell, r)
169						}
170					}
171				} else {
172					if (small) {
173						ppp <- pp@polygons[[1]]@Polygons
174						ishole <- sapply(ppp, function(z)z@hole)
175						xy <- lapply(ppp, function(z)z@coords)
176						xy <- xy[!ishole]
177						if (length(xy) > 0) {
178							cell <- unique(unlist(lapply(xy, function(z) cellFromXY(x, z)), use.names = FALSE))
179							value <- .cellValues(x, cell, layer=layer, nl=nl)
180							if (weights | exact) {
181								weight=rep(1/NROW(value), NROW(value))
182								if (cellnumbers) {
183									r <- cbind(cell, value, weight)
184								} else {
185									r <- cbind(value, weight)
186								}
187							} else if (cellnumbers) {
188								r <- cbind(cell, value)
189							} else {
190								r <- value
191							}
192						} else {
193							r <- NULL
194						}
195					} else {
196						r <- NULL
197					}
198				}
199			}
200			r
201		}
202
203        for (ni in 1:nodes) {
204			.sendCall(cl[[ni]], clFun, list(ni, y[ni,]), tag=ni)
205		}
206
207		for (i in 1:npol) {
208			d <- .recvOneData(cl)
209			if (! d$value$success) {
210				stop('cluster error at polygon: ', i)
211			}
212
213			if (doFun) {
214				if (!is.null(d$value$value)) {
215					if (nl > 1 & !(weights | exact)) {
216						res[[d$value$tag]] <- apply(d$value$value, 2, fun, na.rm=na.rm)
217					} else {
218						res[[d$value$tag]] <- fun(d$value$value, na.rm=na.rm)
219					}
220				}
221			} else {
222				res[[d$value$tag]] <- d$value$value
223			}
224			ni <- ni + 1
225			if (ni <= npol) {
226				.sendCall(cl[[d$node]], clFun, list(ni, y[ni,]), tag=ni)
227			}
228			pbStep(pb, i)
229		}
230
231	} else {
232		for (i in 1:npol) {
233			pp <- y[i,]
234			spbb <- sp::bbox(pp)
235
236			if (spbb[1,1] >= rsbb[1,2] | spbb[1,2] <= rsbb[1,1] | spbb[2,1] >= rsbb[2,2] | spbb[2,2] <= rsbb[2,1]) {
237				# do nothing; res[[i]] <- NULL
238			} else {
239				rc <- crop(rr, extent(pp)+addres)
240				if (exact) {
241					erc <- crop(x, rc)
242					xy <- exactextractr::exact_extract(erc, pp, include_cell=cellnumbers, progress=FALSE)[[1]]
243				} else if (weights) {
244					rc <- .polygonsToRaster(pp, rc, getCover=TRUE, silent=TRUE)
245					rc[rc==0] <- NA
246					xy <- rasterToPoints(rc)
247					if (normalizeWeights) {
248						weight <- xy[,3] / sum(xy[,3])
249					} else {
250						weight <- xy[,3] #/ 100
251					}
252					xy <- xy[,-3,drop=FALSE]
253				} else {
254					rc <- .polygonsToRaster(pp, rc, silent=TRUE)
255					xy <- rasterToPoints(rc)[,-3,drop=FALSE]
256				}
257
258				if (length(xy) > 0)  {  # catch holes or very small polygons
259					if (exact) {
260						if (weights) {
261							if (normalizeWeights) {
262								xy$coverage_fraction <- xy$coverage_fraction / sum(xy$coverage_fraction)
263							}
264							colnames(xy)[ncol(xy)] <- "weight"
265						} else {
266							xy$coverage_fraction  <- NULL
267						}
268						if (cellnumbers) {
269							nms <- colnames(xy)
270							# not good if there is a layer called cell
271							nms <- c("cell", nms[nms != "cell"])
272							xy <- xy[,nms]
273						}
274						if (ncol(xy) == 1) {
275							res[[i]] <- unlist(xy, use.names =FALSE)
276						} else {
277							res[[i]] <- as.matrix(xy)
278						}
279					} else if (weights) {
280						value <- .xyValues(x, xy, layer=layer, nl=nl)
281						if (cellnumbers) {
282							cell <- cellFromXY(x, xy)
283							res[[i]] <- cbind(cell, value, weight)
284						} else {
285							res[[i]] <- cbind(value, weight)
286						}
287					} else if (cellnumbers) {
288						value <- .xyValues(x, xy, layer=layer, nl=nl)
289						cell <- cellFromXY(x, xy)
290						res[[i]] <- cbind(cell, value)
291					} else {
292						res[[i]] <- .xyValues(x, xy, layer=layer, nl=nl)
293					}
294				} else if (small) {
295					ppp <- pp@polygons[[1]]@Polygons
296					ishole <- sapply(ppp, function(z)z@hole)
297					xy <- lapply(ppp, function(z)z@coords)
298					xy <- xy[!ishole]
299					if (length(xy) > 0) {
300						cell <- unique(unlist(lapply(xy, function(z) cellFromXY(x, z))), use.names = FALSE)
301						value <- .cellValues(x, cell, layer=layer, nl=nl)
302						if (weights | exact) {
303							weight <- rep(1/NROW(value), NROW(value))
304							if (cellnumbers) {
305								res[[i]] <- cbind(cell, value, weight)
306							} else {
307								res[[i]] <- cbind(value, weight)
308							}
309						} else if (cellnumbers) {
310							res[[i]] <- cbind(cell, value)
311						} else {
312							res[[i]] <- value
313						}
314					} # else do nothing; res[[i]] <- NULL
315				}
316				if (doFun) {
317					if (!is.null(res[[i]])) {
318						if (nl > 1 & !(weights | exact)) {
319							res[[i]] <- apply(res[[i]], 2, fun, na.rm=na.rm)
320						} else {
321							res[[i]] <- fun(res[[i]], na.rm=na.rm)
322						}
323					}
324				}
325			}
326			pbStep(pb)
327		}
328	}
329	res <- res[1:npol]
330	pbClose(pb)
331
332
333	if (! is.null(fun)) {
334		# try to simplify
335		i <- sapply(res, length)
336		if (length(unique(i[i != 0])) == 1) {
337			if (any(i == 0)) {
338				lng <- length(res)
339				v <- do.call(rbind, res)
340				res <- matrix(NA, nrow=lng, ncol=ncol(v))
341				res[which(i > 0), ] <- v
342			} else {
343				res <- do.call(rbind, res)
344			}
345		} else {
346			if (sp) {
347				warning('cannot return a sp object because the data length varies between polygons')
348				sp <- FALSE
349				df <- FALSE
350			#} else if (df) {
351				#warning('cannot return a data.frame because the data length varies between polygons')
352				#df <- FALSE
353			}
354		}
355	}
356
357	if (df) {
358		if (!is.list(res)) {
359			res <- data.frame(ID=1:NROW(res), res)
360		} else {
361			res <- data.frame( do.call(rbind, lapply(1:length(res), function(x) if (!is.null(res[[x]])) cbind(x, res[[x]]))) )
362		}
363
364		lyrs <- layer:(layer+nl-1)
365		if (cellnumbers) {
366			nms <- c('ID', 'cell', names(x)[lyrs])
367		} else {
368			nms <- c('ID', names(x)[lyrs])
369		}
370		if ((weights|exact) & is.null(fun)) {
371			nms <- c(nms, 'weight')
372		}
373		colnames(res) <- nms
374
375		if (any(is.factor(x)) & factors) {
376			i <- ifelse(cellnumbers, 1:2, 1)
377			v <- res[, -i, drop=FALSE]
378			if (ncol(v) == 1) {
379				v <- data.frame(factorValues(x, v[,1], layer))
380			} else {
381				v <- .insertFacts(x, v, lyrs)
382			}
383			res <- data.frame(res[,i,drop=FALSE], v)
384		}
385	}
386
387	if (sp) {
388		if (nrow(res) != npol) {
389			warning('sp=TRUE is ignored because fun does not summarize the values of each polygon to a single number')
390			return(res)
391		}
392
393		if (!.hasSlot(y, 'data') ) {
394			y <- sp::SpatialPolygonsDataFrame(y, res[, -1, drop=FALSE])
395		} else {
396			y@data <- cbind(y@data, res[, -1, drop=FALSE])
397		}
398		return(y)
399	}
400
401	res
402}
403)
404
405