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