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