1# Author: Robert J. Hijmans 2# Date : March 2009 3# Licence GPL v3 4# updated November 2011 5# version 1.0 6 7 8 9.bilinearValue <- function(raster, xyCoords, layer, n) { 10 11 #bilinear <- function(xy, x, y, v) { 12 # .doBilinear(xy, x, y, v) 13 #} 14 15 r <- raster(raster) 16 nls <- nlayers(raster) 17 18 four <- fourCellsFromXY(r, xyCoords, duplicates=FALSE) 19 20 xy4 <- matrix(xyFromCell(r, as.vector(four)), ncol=8) 21 x <- rbind(.doSpmin(xy4[,1], xy4[,3]), .doSpmax(xy4[,1], xy4[,3])) 22 y <- rbind(.doSpmin(xy4[,5], xy4[,6]), .doSpmax(xy4[,5], xy4[,6])) 23 # data.frame is faster than cbind in this case (less copying?) 24 xy4 <- data.frame( 25 x = c(x[1,], x[1,], x[2,], x[2,]), 26 y = c(y[1,], y[2,], y[1,], y[2,]) 27 ) 28 cells <- cellFromXY(r, xy4) 29 30 suppressWarnings(row1 <- rowFromCell(r, min(cells, na.rm=TRUE))) 31 if (is.na(row1)) { 32 if (nls == 1) { 33 return(rep(NA, nrow(xyCoords))) 34 } else { 35 return(matrix(NA, nrow= nrow(xyCoords), ncol=nls)) 36 } 37 } 38 39 nrows <- rowFromCell(r, max(cells, na.rm=TRUE)) - row1 + 1 40 offs <- cellFromRowCol(r, row1, 1) - 1 41 cells <- cells - offs 42 43 if (nls == 1) { 44 vv <- getValues(raster, row1, nrows) 45 v <- matrix( vv[cells], ncol=4) 46 47 res <- rep(NA, nrow(v)) 48 rs <- rowSums(is.na(v)) 49 i <- rs==3 50 if (sum(i) > 0) { 51 cells <- cellFromXY(raster, xyCoords[i,]) - offs 52 res[i] <- vv[cells] 53 } 54 i <- rs > 0 & rs < 3 55 if (sum(i) > 0) { 56 vv <- v[i,,drop=FALSE] 57 vv[is.na(vv[,1]),1] <- vv[is.na(vv[,1]),2] 58 vv[is.na(vv[,2]),2] <- vv[is.na(vv[,2]),1] 59 vv[is.na(vv[,3]),3] <- vv[is.na(vv[,3]),4] 60 vv[is.na(vv[,4]),4] <- vv[is.na(vv[,4]),3] 61 vmean <- rep(rowMeans(vv, na.rm=TRUE), 4) 62 vv[is.na(vv)] <- vmean[is.na(vv)] 63# res[i] <- bilinear(xyCoords[i,1], xyCoords[i,2], x[1,i], x[2,i], y[1,i], y[2,i], vv) 64 res[i] <- .doBilinear(xyCoords[i,,drop=FALSE], x[,i,drop=FALSE], y[,i,drop=FALSE], vv) 65 } 66 i <- rs==0 67 if (sum(i) > 0) { 68# res[i] <- bilinear(xyCoords[i,1], xyCoords[i,2], x[1,i], x[2,i], y[1,i], y[2,i], v[i,]) 69 res[i] <- .doBilinear(xyCoords[i, ,drop=FALSE], x[,i,drop=FALSE], y[,i,drop=FALSE], v[i,,drop=FALSE]) 70 } 71 res 72 73 } else { 74 75 if (missing(layer)) { layer <- 1 } 76 if (missing(n)) { n <- (nls-layer+1) } 77 lyrs <- layer:(layer+n-1) 78 allres <- matrix(ncol=length(lyrs), nrow=nrow(xyCoords)) 79 colnames(allres) <- names(raster)[lyrs] 80 81 cvv <- getValues(raster, row1, nrows)[, lyrs] 82 cv <- cvv[cells,] 83 for (j in 1:ncol(cv)) { 84 v <- matrix(cv[, j], ncol=4) 85 86 res <- rep(NA, nrow(v)) 87 rs <- rowSums(is.na(v)) 88 i <- rs==3 89 if (sum(i) > 0) { 90 cells <- cellFromXY(raster, xyCoords[i,]) - offs 91 res[i] <- cvv[cells, j] 92 } 93 i <- rs > 0 & rs < 3 94 if (sum(i) > 0) { 95 vv <- v[i,,drop=FALSE] 96 vv[is.na(vv[,1]),1] <- vv[is.na(vv[,1]),2] 97 vv[is.na(vv[,2]),2] <- vv[is.na(vv[,2]),1] 98 vv[is.na(vv[,3]),3] <- vv[is.na(vv[,3]),4] 99 vv[is.na(vv[,4]),4] <- vv[is.na(vv[,4]),3] 100 vmean <- rep(rowMeans(vv, na.rm=TRUE), 4) 101 vv[is.na(vv)] <- vmean[is.na(vv)] 102 res[i] <- .doBilinear(xyCoords[i,,drop=FALSE], x[,i,drop=FALSE], y[,i,drop=FALSE], vv) 103 } 104 i <- rs==0 105 if (sum(i) > 0) { 106 res[i] <- .doBilinear(xyCoords[i,,drop=FALSE], x[,i,drop=FALSE], y[,i,drop=FALSE], v[i,,drop=FALSE]) 107 } 108 109 allres[,j] <- res 110 } 111 allres 112 } 113} 114 115