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