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