1# Author: Robert J. Hijmans
2# Date :  November 2009
3# Version 1.0
4# Licence GPL v3
5
6# under development
7
8..idwValue <- function(raster, xy, ngb=4, pow=1, layer, n) {
9	r <- raster(raster)
10	longlat <- couldBeLonLat(r)
11	cells <- cellFromXY(r, xy)
12	adj <- adjacent(r, cells, ngb, pairs=TRUE, include=TRUE, id=TRUE)
13
14	uc <- unique(adj[,3])
15	row1 <- rowFromCell(r, min(uc, na.rm=TRUE))
16	nrows <- row1 - 1 + rowFromCell(r, max(uc, na.rm=TRUE))
17	offs <- cellFromRowCol(r, row1, 1) - 1
18	cs <- uc - offs
19
20	nl <- nlayers(raster)
21	if (nl==1) {
22		v <- cbind(uc, v=getValues(raster, row1, nrows)[cs])
23	} else {
24		v <- cbind(uc, v=getValues(raster, row1, nrows)[cs,])
25	}
26	m <- merge(adj, v, by.x='to', by.y=1)
27	colnames(xy) <- c('x', 'y')
28	m <- merge(m, cbind(1:nrow(xy), xy), by.x='id', by.y=1)
29
30	pd <- pointDistance(m[,c('x', 'y')], xyFromCell(r, m$to), lonlat=longlat) / 1000
31	pd <- pd^pow
32	pd[pd==0] <- 1e-12
33
34	if (nl==1) {
35		pd[is.na(m$v)] <- NA
36		as.vector( tapply(m$v*(1/pd), m$id, sum, na.rm=TRUE) / tapply(1/pd, m$id, sum, na.rm=TRUE) )
37		#cbind(as.integer(names(res)), res)
38	} else {
39		lys <- 4:(4+nl-1)
40		a1 <- aggregate(m[,lys]*(1/pd), list(m$id), sum)
41		a2 <- aggregate(1/pd, list(m$id), sum)
42		res <- as.matrix(a1[,-1]) / as.vector(as.matrix(a2[,-1]))
43		res <- cbind(as.vector(a1[,1]), res)
44		res[, -1]
45	}
46}
47
48#a=raster(nc=10,nr=10)
49#xmin(a)=55
50#projection(a) = "+proj=utm +zone=33"
51#a[] = 1:ncell(a)
52#a[50:75]=NA
53#r = disaggregate(raster(a), 3)
54#r[] = .idwValue(a, sp::coordinates(r))
55#plot(r)
56