1# Author: Robert J. Hijmans 2# Date : January 2009 3# Version 2.0 4# Licence GPL v3 5 6 7.getPutVals <- function(obj, field, n, mask) { 8 9 if (mask) { 10 return( data.frame(v=rep(1, length=n)) ) 11 12 } else if (missing(field)) { 13 if (.hasSlot(obj, 'data')) { 14 putvals <- obj@data 15 cn <- validNames(c('ID', colnames(putvals))) 16 cn[1] <- 'ID' 17 putvals <- data.frame(ID=1:nrow(putvals), putvals) 18 colnames(putvals) <- cn 19 } else { 20 putvals <- data.frame(v=as.integer(1:n)) 21 } 22 return(putvals) 23 24 } else if (isTRUE (is.na(field))) { 25 return( data.frame(v=rep(NA, n)) ) 26 27 28 } else if (is.character(field) ) { 29 if (.hasSlot(obj, 'data')) { 30 nms <- names(obj) 31 if (length(field) <= length(nms)) { 32 m <- match(field, nms) 33 if (!all(is.na(m))) { 34 m <- stats::na.omit(m) 35 return(obj@data[, m, drop=FALSE]) 36 } 37 } 38 } 39 } 40 41 if (NROW(field) == n) { 42 if (is.null(nrow(field))) { 43 return(data.frame(field, stringsAsFactors=FALSE)) 44 } else { 45 return(field) 46 } 47 48 } 49 50 if (is.numeric(field)) { 51 putvals <- rep(field, length.out=n) 52 return(data.frame(field=putvals)) 53 } 54 55 stop('invalid value for field') 56} 57 58 59.intersectSegments <- function(x1, y1, x2, y2, x3, y3, x4, y4) { 60# Translated by RH from LISP code by Paul Reiners 61# http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/linesegments.lisp 62# Which was translated from the algorithm by Paul Bourke given here: http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/ 63 denom <- ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1)) 64 ua_num <- ((x4 - x3) *(y1 - y3)) - ((y4 - y3) * (x1 - x3)) 65 ub_num <- ((x2 - x1) *(y1 - y3)) - ((y2 - y1) * (x1 - x3)) 66# If the denominator and numerator for the equations for ua and ub are 0 then the two lines are coincident. 67 if ( denom == 0 ) { 68 if (ua_num == 0 & ub_num == 0) { 69 xmin <- max(x1, x3) 70 if (xmin==x1) {ymin <- y1} else {ymin <- y3} 71 xmax <- min(x2, x4) 72 if (xmax==x2) {ymax <- y2} else {ymax <- y4} 73 # RH: for coincident line (segments) returning two intersections : start and end 74 return(rbind(c(xmin, ymin), 75 c(xmax, ymax))) 76 } #else { 77# If the denominator for the equations for ua and ub is 0 then the two lines are parallel. 78# return(c(NA, NA)) 79# } 80 } else { 81 ua <- round(ua_num / denom, 12) 82 ub <- round(ub_num / denom, 12) 83 if ((ua >= 0 & ua <= 1) & (ub >= 0 & ub <= 1) ) { 84 x <- x1 + ua * (x2 - x1) 85 y <- y1 + ua * (y2 - y1) 86 return(c(x, y)) 87 } 88 } 89 90 return(c(NA, NA)) 91} 92 93 94.intersectLinePolygon <- function(line, poly) { 95 resxy <- matrix(NA, ncol=2, nrow=0) 96 miny <- min(line[,2]) 97 maxy <- max(line[,2]) 98 xyxy <- cbind(poly, rbind(poly[-1,], poly[1,])) 99 xyxy <- subset(xyxy, !( (xyxy[,2] > maxy & xyxy[,4] > maxy ) | (xyxy[,2] < miny & xyxy[,4] < miny)) ) 100 if (nrow(xyxy) == 0) { 101 return(resxy) 102 } 103 for (i in 1:nrow(xyxy)) { 104 xy <- .intersectSegments(xyxy[i,1], xyxy[i,2], xyxy[i,3], xyxy[i,4], line[1,1], line[1,2], line[2,1], line[2,2] ) 105 if (!is.na(xy[1])) { 106 resxy <- rbind(resxy, xy) 107 } 108 } 109 return((resxy)) 110} 111 112 113 114 115 116.polygonsToRaster <- function(p, rstr, field, fun='last', background=NA, mask=FALSE, update=FALSE, updateValue="all", getCover=FALSE, filename="", silent=TRUE, faster=TRUE, ...) { 117 118 119 npol <- length(p@polygons) 120 pvals <- .getPutVals(p, field, npol, mask) 121 putvals <- pvals[,1] 122 if (ncol(pvals) > 1) { 123 rstr@data@isfactor <- TRUE 124 rstr@data@attributes <- list(pvals) 125 if (!is.character(fun)) { 126 stop('when rasterizing multiple fields you must use "fun=first" or "fun=last"') 127 } else if (!(fun %in% c('first', 'last'))) { 128 stop('when rasterizing multiple fields you must use "fun=first" or "fun=last"') 129 } 130 } 131 132 133 if (getCover) { 134 nc <- ncell(rstr) 135 # high precision for possibly small polygons 136 #https://stackoverflow.com/questions/53854910/issue-with-estimating-weighted-mean-from-raster-for-a-polygon-shape-in-r/ 137 fctr <- ifelse(nc < 5, 100, ifelse(nc < 17, 20, 10)) 138 rstr <- disaggregate(raster(rstr), fctr) 139 r <- .fasterize(p, rstr, rep(1, npol), background=0, datatype="INT1U") 140 return( aggregate(r, fctr, mean, na.rm=TRUE, filename=filename, ...) ) 141 } 142 143 144 145 ### new code 146 if (is.character(fun) && (ncol(pvals) == 1) && faster) { 147 148 if (fun == "last") { 149 if (mask || update) { 150 if (mask && update) stop("either use 'mask' OR 'update'") 151 background = NA 152 r <- .fasterize(p, rstr, pvals[,1], background) 153 if (! hasValues(r)) { 154 if (mask) { 155 warning('there are no values to mask') 156 } else { 157 warning('there are no values to update') 158 } 159 return(r) 160 } 161 if (mask) { 162 r <- mask(rstr, r) 163 } else { 164 if (updateValue[1]=="all") { 165 r <- cover(r, rstr) 166 } else if (updateValue[1]=="NA") { 167 r <- cover(rstr, r, ...) 168 } else if (updateValue[1]=="!NA") { 169 r <- mask(cover(r, rstr), rstr, ...) 170 } else { 171 s <- stack(r, rstr) 172 r <- overlay(rstr, r, fun=function(x,y){ i = (x %in% updateValue & !is.na(y)); x[i] <- y[i]; x }, ... ) 173 } 174 } 175 return(r) 176 } else { 177 return( .fasterize(p, rstr, pvals[,1], background, filename, ...) ) 178 } 179 } 180 181 } 182 ### end new code 183 184 185 186 leftColFromX <- function ( object, x ) { 187 colnr <- (x - xmin(object)) / xres(object) 188 i <- colnr %% 1 == 0 189 colnr[!i] <- trunc(colnr[!i]) + 1 190 colnr[colnr <= 0] <- 1 191 colnr 192 } 193 194 195 rightColFromX <- function ( object, x ) { 196 colnr <- trunc((x - xmin(object)) / xres(object)) + 1 197 colnr[ colnr > ncol(object) ] <- object@ncols 198 colnr 199 } 200 201 202 if (! inherits(p, 'SpatialPolygons') ) { 203 stop('The first argument should be an object of the "SpatialPolygons*" lineage') 204 } 205 206 filename <- trim(filename) 207 if (!canProcessInMemory(rstr, 3) && filename == '') { 208 filename <- rasterTmpFile() 209 } 210 211 212 213 if (mask & update) { 214 stop('use either "mask" OR "update"') 215 } else if (mask) { 216 oldraster <- rstr 217 #update <- TRUE 218 } else if (update) { 219 oldraster <- rstr 220 if (!is.numeric(updateValue)) { 221 if (is.na(updateValue)) { 222 updateValue <- 'NA' 223 } else if (!(updateValue == 'NA' | updateValue == '!NA' | updateValue == 'all')) { 224 stop('updateValue should be either "all", "NA", "!NA"') 225 } 226 } 227 } 228 229 rstr <- raster(rstr) 230 231 if (!is.na(projection(p))) { 232 projection(rstr) <-.getCRS(p) 233 } 234 235# check if bbox of raster and p overlap 236 spbb <- sp::bbox(p) 237 rsbb <- bbox(rstr) 238 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]) { 239 # instead of a warning 240 return( init(rstr, function(x) NA) ) 241 # so that clusterR can use this function (overlap with some chunks might be NULL) 242 } 243 244 npol <- length(p@polygons) 245 pvals <- .getPutVals(p, field, npol, mask) 246 putvals <- pvals[,1] 247 if (ncol(pvals) > 1) { 248 rstr@data@isfactor <- TRUE 249 rstr@data@attributes <- list(pvals) 250 if (!is.character(fun)) { 251 stop('when rasterizing multiple values you must use "fun=first" or "fun=last"') 252 } else if (!(fun %in% c('first', 'last'))) { 253 stop('when rasterizing multiple values you must use "fun=first" or "fun=last"') 254 } 255 } 256 257 if (is.character(fun)) { 258 if (fun=='first') { 259 fun <- function(x, ...){ stats::na.omit(x)[1] } 260 } else if (fun=='last') { 261 fun <- function(x, ...){ rev(stats::na.omit(x))[1] } 262 } else if (fun == 'count') { 263 fun <- function(x, ...){ sum(!is.na(x)) } 264 field <- 1 265 } 266 } 267 268 polinfo <- data.frame(matrix(NA, nrow=npol * 2, ncol=6)) 269 colnames(polinfo) <- c('part', 'miny', 'maxy', 'value', 'hole', 'object') 270 addpol <- polinfo[rep(1, 500), ] 271 rownames(addpol) <- NULL 272 pollist <- list() 273 cnt <- 0 274 for (i in 1:npol) { 275 nsubpol <- length(p@polygons[[i]]@Polygons) 276 for (j in 1:nsubpol) { 277 cnt <- cnt + 1 278 if (cnt > dim(polinfo)[1]) { 279 polinfo <- rbind(polinfo, addpol) 280 } 281 polinfo[cnt, 1] <- cnt 282 polinfo[cnt, 2] <- min(p@polygons[[i]]@Polygons[[j]]@coords[,2]) 283 polinfo[cnt, 3] <- max(p@polygons[[i]]@Polygons[[j]]@coords[,2]) 284 polinfo[cnt, 4] <- putvals[i] 285 if ( p@polygons[[i]]@Polygons[[j]]@hole ) { 286 polinfo[cnt, 5] <- 1 287 } else { 288 polinfo[cnt, 5] <- 0 289 } 290 polinfo[cnt, 6] <- i 291 pollist[[cnt]] <- p@polygons[[i]]@Polygons[[j]] 292 } 293 } 294 295 if (! silent) { 296 message('Found ', npol, ' region(s) and ', cnt, ' polygon(s)') 297 } 298 299 polinfo <- subset(polinfo, polinfo[,1] <= cnt, drop=FALSE) 300# polinfo <- polinfo[order(polinfo[,1]),] 301# rm(p) 302 303 lxmin <- min(spbb[1,1], rsbb[1,1]) - xres(rstr) 304 lxmax <- max(spbb[1,2], rsbb[1,2]) + xres(rstr) 305 306# if (getCover) { 307# return (.polygoncover(rstr, filename, polinfo, lxmin, lxmax, pollist, ...)) 308# } 309 310 adj <- 0.5 * xres(rstr) 311 312 if (filename == "") { 313 v <- matrix(NA, ncol=nrow(rstr), nrow=ncol(rstr)) 314 } else { 315 rstr <- writeStart(rstr, filename=filename, ...) 316 } 317 318 rxmn <- xmin(rstr) 319 rxmx <- xmax(rstr) 320 321 rv1 <- rep(NA, ncol(rstr)) 322 holes1 <- rep(0, ncol(rstr)) 323 324 pb <- pbCreate(nrow(rstr), label='rasterize', ...) 325 326 for (r in 1:nrow(rstr)) { 327 328 vals <- NULL 329 holes <- holes1 330 331 ly <- yFromRow(rstr, r) 332 myline <- rbind(c(lxmin,ly), c(lxmax,ly)) 333 334 subpol <- subset(polinfo, !(polinfo[,2] > ly | polinfo[,3] < ly), drop=FALSE) 335 if (length(subpol[,1]) > 0) { 336 updateHoles <- FALSE 337 lastpolnr <- subpol[1,6] 338 rvtmp <- rv1 339 for (i in 1:nrow(subpol)) { 340 if (i == nrow(subpol)) { 341 updateHoles <- TRUE 342 } else if (subpol[i+1,6] > lastpolnr) { # new polygon 343 updateHoles <- TRUE 344 lastpolnr <- subpol[i+1,6] 345 } 346 347 mypoly <- pollist[[subpol[i,1]]] 348 intersection <- .intersectLinePolygon(myline, mypoly@coords) 349 #if (nrow(intersection) %% 2 == 1) { 350 # this is a bit speculative 351 # not OK! 352 # intersection <- unique(intersection) 353 #} 354 355 x <- sort(intersection[,1]) 356 if (length(x) > 0) { 357 if ((nrow(intersection) %% 2 == 1) || ( sum(x[-length(x)] == x[-1]) > 0 )) { 358 # uneven number or duplicates 359 # e.g. single node intersection going out of polygon .... 360 spPnts <- sp::SpatialPoints(xyFromCell(rstr, cellFromRowCol(rstr, rep(r, ncol(rstr)), 1:ncol(rstr)))) 361 spPol <- sp::SpatialPolygons(list(sp::Polygons(list(mypoly), 1))) 362 over <- sp::over(spPnts, spPol) 363 if ( subpol[i, 5] == 1 ) { 364 holes[!is.na(over)] <- holes[!is.na(over)] - 1 365 } else { 366 rvtmp[!is.na(over)] <- subpol[i,4] 367 holes[!is.na(over)] <- holes[!is.na(over)] + 1 368 } 369 # print(paste('exit node intersection on row:', r)) 370 } else { 371 372 for (k in 1:round(nrow(intersection)/2)) { 373 l <- (k * 2) - 1 374 x1 <- x[l] 375 x2 <- x[l+1] 376 #if (is.na(x2)) { 377 # txt <- paste('something funny at row:', r, 'polygon:',j) 378 # stop(txt) 379 #} 380 # if (x1 > rxmx) { next } 381 # if (x2 < rxmn) { next } 382 # adjust to skip first cell if the center is not covered by this polygon 383 x1a <- x1 + adj 384 x2a <- x2 - adj 385 if (x1a > rxmx) { next } 386 if (x2a < rxmn) { next } 387 x1a <- min(rxmx, max(rxmn, x1a)) 388 x2a <- min(rxmx, max(rxmn, x2a)) 389 col1 <- leftColFromX(rstr, x1a) 390 col2 <- rightColFromX(rstr, x2a) 391 if (col1 > col2) { 392 spPnts <- sp::SpatialPoints(xyFromCell(rstr, cellFromRowCol(rstr, rep(r, ncol(rstr)), 1:ncol(rstr)))) 393 spPol <- sp::SpatialPolygons(list(sp::Polygons(list(mypoly), 1))) 394 over <- sp::over(spPnts, spPol) 395 if ( subpol[i, 5] == 1 ) { 396 holes[!is.na(over)] <- holes[!is.na(over)] - 1 397 } else { 398 rvtmp[!is.na(over)] <- subpol[i,4] 399 holes[!is.na(over)] <- holes[!is.na(over)] + 1 400 } 401 next 402 } 403 if ( subpol[i, 5] == 1 ) { 404 holes[col1:col2] <- holes[col1:col2] - 1 405 } else { 406 rvtmp[col1:col2] <- subpol[i,4] 407 holes[col1:col2] <- holes[col1:col2] + 1 408 } 409 } 410 } 411 } 412 413 if (updateHoles) { 414 updateHoles <- FALSE 415 rvtmp[holes < 1] <- NA 416 vals <- cbind(vals, rvtmp) 417 rvtmp <- rv1 418 holes <- holes1 419 } 420 } 421 } 422 423 #print(vals) 424 425 rrv <- rv1 426 if (!is.null(vals)) { 427 u <- which(rowSums(is.na(vals)) < ncol(vals)) 428 if (length(u) > 0) { 429 if (mask) { 430 rrv[u] <- 1 431 } else { 432 rrv[u] <- apply(vals[u, ,drop=FALSE], 1, fun, na.rm=TRUE) 433 } 434 } 435 } 436 437 if (mask) { 438 oldvals <- getValues(oldraster, r) 439 ind <- which(is.na(rrv)) 440 oldvals[ind] <- NA 441 rrv <- oldvals 442 } else if (update) { 443 oldvals <- getValues(oldraster, r) 444 if (is.numeric(updateValue)) { 445 ind <- which(oldvals == updateValue & !is.na(rrv)) 446 } else if (updateValue == "all") { 447 ind <- which(!is.na(rrv)) 448 } else if (updateValue == "NA") { 449 ind <- which(is.na(oldvals)) 450 } else { "!NA" 451 ind <- which(!is.na(oldvals) & !is.na(rrv)) 452 } 453 oldvals[ind] <- rrv[ind] 454 rrv <- oldvals 455 } else { 456 rrv[is.na(rrv)] <- background 457 } 458 459 if (filename == "") { 460 v[,r] <- rrv 461 } else { 462# print(rrv) 463 rstr <- writeValues(rstr, rrv, r) 464 } 465 pbStep(pb, r) 466 } 467 pbClose(pb) 468 469 if (filename == "") { 470 rstr <- setValues(rstr, as.vector(v)) 471 } else { 472 rstr <- writeStop(rstr) 473 } 474 return(rstr) 475} 476 477#plot( .polygonsToRaster(p, rstr) ) 478 479 480#...polygoncover <- function(p, x, filename, ...) { 481# d <- disaggregate(raster(x), 10) 482# r <- .polygonsToRaster(p, d, filename=filename, field=1, fun='first', background=0, mask=FALSE, update=FALSE, getCover=FALSE, silent=TRUE, ...) 483# aggregate(r, 10, sum) 484#} 485 486.Old_polygoncover <- function(rstr, filename, polinfo, lxmin, lxmax, pollist, ...) { 487# percentage cover per grid cell 488 489 polinfo[, 4] <- 1 490 491 bigraster <- raster(rstr) 492 rxmn <- xmin(bigraster) 493 rxmx <- xmax(bigraster) 494 f <- 10 495 adj <- 0.5 * xres(bigraster)/f 496 nc <- ncol(bigraster) * f 497 rv1 <- rep(0, nc) 498 holes1 <- rep(0, nc) 499 prj <-.getCRS(bigraster) 500 hr <- 0.5 * yres(bigraster) 501 502 vv <- matrix(ncol=f, nrow=nc) 503 504 if (filename == "") { 505 v <- matrix(NA, ncol=nrow(bigraster), nrow=ncol(bigraster)) 506 } else { 507 bigraster <- writeStart(bigraster, filename=filename, ...) 508 } 509 510 pb <- pbCreate(nrow(bigraster), label='rasterize', ...) 511 for (rr in 1:nrow(bigraster)) { 512 y <- yFromRow(bigraster, rr) 513 yn <- y - hr 514 yx <- y + hr 515 rstr <- raster(xmn=rxmn, xmx=rxmx, ymn=yn, ymx=yx, ncols=nc, nrows=f, crs=prj) 516 subpol <- subset(polinfo, !(polinfo[,2] > yx | polinfo[,3] < yn), drop=FALSE) 517 for (r in 1:f) { 518 rv <- rv1 519 ly <- yFromRow(rstr, r) 520 myline <- rbind(c(lxmin,ly), c(lxmax,ly)) 521 holes <- holes1 522 if (length(subpol[,1]) > 0) { 523 524 updateHoles <- FALSE 525 lastpolnr <- subpol[1,6] 526 rvtmp <- rv1 527 for (i in 1:length(subpol[,1])) { 528 if (i == length(subpol[,1])) { 529 updateHoles <- TRUE 530 } else if (subpol[i+1,6] > lastpolnr) { 531 updateHoles <- TRUE 532 lastpolnr <- subpol[i+1,6] 533 } 534 535 mypoly <- pollist[[subpol[i,1]]] 536 intersection <- .intersectLinePolygon(myline, mypoly@coords) 537 x <- sort(intersection[,1]) 538 if (length(x) > 0) { 539 #if (length(subpol[,1]) > 3 & i ==2) { 540 # print('4') 541 #} 542 if ( sum(x[-length(x)] == x[-1]) > 0 ) { 543 # single node intersection going out of polygon .... 544 spPnts <- sp::SpatialPoints(xyFromCell(rstr, cellFromRowCol(rstr, rep(r, ncol(rstr)), 1:ncol(rstr)))) 545 spPol <- sp::SpatialPolygons(list(sp::Polygons(list(mypoly), 1))) 546 over <- sp::over(spPnts, spPol) 547 if ( subpol[i, 5] == 1 ) { 548 holes[!is.na(over)] <- holes[!is.na(over)] - 1 549 } else { 550 rvtmp[!is.na(over)] <- subpol[i,4] 551 holes[!is.na(over)] <- holes[!is.na(over)] + 1 552 } 553 } else { 554 for (k in 1:round(nrow(intersection)/2)) { 555 l <- (k * 2) - 1 556 x1 <- x[l] 557 x2 <- x[l+1] 558 if (x1 > rxmx) { next } 559 if (x2 < rxmn) { next } 560 # adjust to skip first cell if the center is not covered by this polygon 561 x1a <- x1 + adj 562 x2a <- x2 - adj 563 x1a <- min(rxmx, max(rxmn, x1a)) 564 x2a <- min(rxmx, max(rxmn, x2a)) 565 col1 <- colFromX(rstr, x1a) 566 col2 <- colFromX(rstr, x2a) 567 if (col1 > col2) { next } 568 if ( subpol[i, 5] == 1 ) { 569 holes[col1:col2] <- holes[col1:col2] - 1 570 } else { 571 rvtmp[col1:col2] <- subpol[i,4] 572 holes[col1:col2] <- holes[col1:col2] + 1 573 } 574 } 575 } 576 if (updateHoles) { 577 holes <- holes < 1 578 rvtmp[holes] <- 0 579 holes <- holes1 580 updateHoles <- FALSE 581 rv <- pmax(rv, rvtmp) 582 } 583 } 584 } 585 } 586 vv[,r] <- rv 587 } 588 av <- colSums( matrix( rowSums(vv), nrow=f) ) 589 590 if (filename == "") { 591 v[,rr] <- av 592 } else { 593 bigraster <- writeValues(bigraster, av, rr) 594 } 595 pbStep(pb, rr) 596 } 597 pbClose(pb) 598 599 if (filename == "") { 600 bigraster <- setValues(bigraster, as.vector(v)) 601 } else { 602 bigraster <- writeStop(bigraster) 603 } 604 return(bigraster) 605} 606 607 608#x = .polygoncover(rstr, "", polinfo, lxmin, lxmax, pollist) 609 610 611.polygonsToRaster2 <- function(p, raster, field=0, filename="", ...) { 612# This is based on sampling by points. Should be slower except when polygons very detailed and raster has low resolution 613# but it could be optimized further 614# currently not used. Perhaps it should be used under certain conditions. 615# this version does not deal with polygon holes 616 617# check if bbox of raster and p overlap 618 filename <- trim(filename) 619 raster <- raster(raster) 620 621 spbb <- sp::bbox(p) 622 rsbb <- bbox(raster) 623 if (spbb[1,1] > rsbb[1,2] | spbb[2,1] > rsbb[2,2]) { 624 stop('polygon and raster have no overlapping areas') 625 } 626 627 if (class(p) == 'SpatialPolygons' | field == 0) { 628 putvals <- 1:length(p@polygons) 629 } else { 630 putvals <- as.vector(p@data[,field]) 631 if (class(putvals) == 'character') { 632 stop('selected field is charater type') 633 } 634 } 635 636 637 if (filename == "") { 638 v <- vector(length=0) # replace this 639 } else { 640 raster <- writeStart(raster, filename=filename, ...) 641 } 642 643 rowcol <- cbind(0, 1:ncol(raster)) 644 645 firstrow <- rowFromY(raster, spbb[2,2]) 646 lastrow <- rowFromY(raster, spbb[2,1]) 647 648 for (r in 1:nrow(raster)) { 649 if (r < firstrow | r > lastrow) { 650 vals <- rep(NA, times=ncol(raster)) 651 } else { 652 rowcol[,1] <- r 653 sppoints <- xyFromCell(raster, cellFromRowCol(raster, rowcol[,1], rowcol[,2]), TRUE) 654 over <- sp::over(sppoints, p) 655 vals <- putvals[over] 656 } 657 if (filename == "") { 658 v <- c(v, vals) 659 } else { 660 raster <- writeValues(raster, vals) 661 } 662 } 663 if (filename == "") { 664 raster <- setValues(raster, v) 665 } else { 666 raster <- writeStop(raster) 667 } 668 return(raster) 669} 670 671