# R/blend.R In raster: Geographic Data Analysis and Modeling

#### Defines functions blendold_blend

# Authors: Rafael Wueest, WSL Birmensdorf, Switzerland, [email protected],
# Etienne B. Racine, Robert J. Hijmans
# Date : November 2012
# Version 1.0
# Licence GPL v3

# needs to be generalized to n input rasters and to multi-layer objects
.old_blend <- function(r1, r2) {
i <- intersect(raster(r1), raster(r2))
j <- extend(i, c(1,1))
a <- crop(r1, j)
b <- crop(r2, j)
values(a) <- 1
values(b) <- 2
ab <- merge(a, b)
ba <- merge(b, a)
p1 <- rasterToPoints(ab, function(x) x==2)
p2 <- rasterToPoints(ba, function(x) x==1)
d1 <- distanceFromPoints(i, p1[,1:2])
d2 <- distanceFromPoints(i, p2[,1:2])
dsum <- d1 + d2

z1 <- d1 * crop(r1, d1) / dsum
z2 <- d2 * crop(r2, d2) / dsum
merge(z1 + z2, r1, r2)
}

.blend <- function(x, y, logistic=FALSE, filename='', ...) {

# check for difference in extent
stopifnot( extent(x) != extent(y))

# define logistic function
if (logistic) {
G <- 1
f <- 0.001
k <- log(G/f-1)/(0.5*G)
logfun <- function(x) {
G /(1+exp(-k*G*x)*(G/f-1))
}
}

# create intersection rasters
i <- intersect(raster(x), raster(y))
j <- extend(i, c(1,1))

# is one of the rasters nested within the other?
ex <- extent(x)
ey <- extent(y)
exy <- union(ex, ey)

if (exy==ex | exy==ey){    # the nested case

# which raster has the smaller extent?
if (extent(x) < extent(y)){
rlarge <- y
rsmall <- x
} else {
rlarge <- x
rsmall <- y
}

# create points around nested raster
a <- crop(rlarge, j)
a <- setValues(a, 1)
b <- crop(rsmall, j)
b <- setValues(b, 2)
ba <- merge(b, a)
p <- rasterToPoints(ba, function(x) x==1)

# calculate distances to points in nested raster
d <- distanceFromPoints(i, p[,1:2])

# standardize these distances
dmin <- cellStats(d,'min')
d.sc <- (d - dmin + 1e-9) / (cellStats(d,'max') - dmin)

# the logistic case
if(logistic){
d.sc<-logfun(d.sc)
}

# create distance weighted rasters
z1 <- d.sc * crop(rsmall, d.sc)
z2 <- abs(1-d.sc) * crop(rlarge, d.sc)

# merge rasters
m <- merge(z1 + z2, rsmall, rlarge, filename=filename, ...)

} else {    # the overlapping case

# create points around ovelapping area
a <- crop(x, j)
a <- setValues(a, 1)
b <- crop(y, j)
b <- setValues(b, 2)
ab <- merge(a, b)
ba <- merge(b, a)
p1 <- rasterToPoints(ab, function(x) x==2)
p2 <- rasterToPoints(ba, function(x) x==1)

# calculate distances to points in overlapping area
d1 <- distanceFromPoints(i, p1[,1:2])
d2 <- distanceFromPoints(i, p2[,1:2])

# the logistic case
if(logistic){
d1min <- cellStats(d1,'min')
d2min <- cellStats(d2,'min')
d1 <- logfun((d1 - d1min + 1e-9)/(cellStats(d1,'max') - d1min))
d2 <- logfun((d2 - d2min + 1e-9)/(cellStats(d2,'max') - d2min))
}

# sum distance rasters
dsum <- d1 + d2

# create distance weighted rasters
z1 <- d1 * crop(x, d1) / dsum
z2 <- d2 * crop(y, d2) / dsum
z <- sum(z1, z2)

# merge rasters
m <- merge(z, x, y, filename=filename, ...)
}

m
}

## Try the raster package in your browser

Any scripts or data that you put into this service are public.

raster documentation built on April 6, 2019, 3:05 a.m.