# Conditional assignment of values to adjacent raster cells?

I have a value raster:

``````m <- matrix(c(2,4,5,5,2,8,7,3,1,6,
5,7,5,7,1,6,7,2,6,3,
4,7,3,4,5,3,7,9,3,8,
9,3,6,8,3,4,7,3,7,8,
3,3,7,7,5,3,2,8,9,8,
7,6,2,6,5,2,2,7,7,7,
4,7,2,5,7,7,7,3,3,5,
7,6,7,5,9,6,5,2,3,2,
4,9,2,5,5,8,3,3,1,2,
5,2,6,5,1,5,3,7,7,2),nrow=10, ncol=10, byrow = T)
r <- raster(m)
extent(r) <- matrix(c(0, 0, 10, 10), nrow=2)
plot(r)
text(r)
``````

From this raster, how can I assign values (or change values) to the 8 adjacent cells of the current cell according to this illustration ? I placed a red point within the current cell from this code line:

``````points(xFromCol(r, col=5), yFromRow(r, row=5),col="red",pch=16)
`````` Here, the expected result will be: where the value of the current cell (i.e, 5 in the value raster) is replaced with 0.

Overall, the new values for the 8 adjacent cells must be calculated as follows:

New value = average of cell values contained in the red rectangle * distance between the current cell (red point) and the adjacent cell (i.e., sqrt(2) for diagonally adjacent cells or 1 otherwise)

Update

When bounds for the adjacent cells are out of the raster limits, I need to calculate new values for the adjacent cells which respect the conditions. The adjacent cells which don't respect the conditions will equal to "NA".

For example, if the reference position is c(1,1) instead of c(5,5) by using [row, col] notation, only the new value at the bottom-right corner can be calculated. Thus, the expected result will be:

``````     [,1] [,2] [,3]
[1,] NA   NA   NA
[2,] NA   0    NA
[3,] NA   NA   New_value
``````

For example, if the reference position is c(3,1), only the new values at the top-right, right and bottom-right corners can be calculated. Thus, the expected result will be:

``````     [,1] [,2] [,3]
[1,] NA   NA   New_value
[2,] NA   0    New_value
[3,] NA   NA   New_value
``````

Here is my first attempt at this by using the function `focal` but I have some difficulty to make an automatic code.

Select adjacent cells

``````mat_perc <- matrix(c(1,1,1,1,1,
1,1,1,1,1,
1,1,0,1,1,
1,1,1,1,1,
1,1,1,1,1), nrow=5, ncol=5, byrow = T)
cell_perc <- adjacent(r, cellFromRowCol(r, 5, 5), directions=mat_perc, pairs=FALSE, sorted=TRUE, include=TRUE)
r_perc <- rasterFromCells(r, cell_perc)
r_perc <- setValues(r_perc,extract(r, cell_perc))
plot(r_perc)
text(r_perc)
``````

if the adjacent cell is located at the upper-left corner of the current cell

``````focal_m <- matrix(c(1,1,NA,1,1,NA,NA,NA,NA), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)*sqrt(2)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the upper-middle corner of the current cell

``````focal_m <- matrix(c(1,1,1,1,1,1,NA,NA,NA), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the upper-left corner of the current cell

``````focal_m <- matrix(c(NA,1,1,NA,1,1,NA,NA,NA), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)*sqrt(2)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the left corner of the current cell

``````focal_m <- matrix(c(1,1,NA,1,1,NA,1,1,NA), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the right corner of the current cell

``````focal_m <- matrix(c(NA,1,1,NA,1,1,NA,1,1), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the bottom-left corner of the current cell

``````focal_m <- matrix(c(NA,NA,NA,1,1,NA,1,1,NA), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)*sqrt(2)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the bottom-middle corner of the current cell

``````focal_m <- matrix(c(NA,NA,NA,1,1,1,1,1,1), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````

if the adjacent cell is located at the bottom-right corner of the current cell

``````focal_m <- matrix(c(NA,NA,NA,NA,1,1,NA,1,1), nrow=3, ncol=3, byrow = T)
focal_function <- function(x) mean(x,na.rm=T)*sqrt(2)
test <- as.matrix(focal(r_perc, focal_m, focal_function))
``````
• +1 I wish all questions were this well framed! Are you looking for a focal operation (moving window statistics)? Check out R's `raster` package and the `focal()` function (p. 90 documentation): cran.r-project.org/web/packages/raster/raster.pdf – Aaron Sep 21 '16 at 23:37
• Thanks very much Aaron for your advice ! Indeed, the function focal seems to be very useful but I am not familiar with it. For example, for the adjacent cell = 8 (figure at the top-left corner), I tested `mat <- matrix(c(1,1,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0), nrow=5, ncol=5, byrow = T)` `f.rast <- function(x) mean(x)*sqrt(2)` `aggr <- as.matrix(focal(r, mat, f.rast))`. How can I obtain the result for only the 8 adjacent cells of the current cell and not all the raster ? Here, the result should be: `res <- matrix(c(7.42,0,0,0,0,0,0,0,0), nrow=3, ncol=3, byrow = T)`. Thanks a lot ! – Pierre Sep 22 '16 at 16:27
• @Pierre Do you need to calculate adjacent values only for position row 5, col 5? Or move this reference position for example to a new reference position row 6, col 6? – Guzmán Sep 26 '16 at 14:30
• Can you explain more (editing your question) about how you need to calculate the adjacent values when the bounds for the adjacent cells are out of the raster limits? E.g.: row 1, col 1. – Guzmán Sep 26 '16 at 17:02
• You examples don't make sense. In the first one, if the reference position is c(1,1), then only the bottom right c(2,2) will get the new value but you have showed that c(3,3) is getting the New_Value. In addition the c(1,1) will become 0 not c(2,2). – Farid Cheraghi Sep 28 '16 at 12:05

## 3 Answers

The function `AssignValuesToAdjacentRasterCells` below returns a new RasterLayer object with the desired values assigned from the original raster input. The function check if the adjacent cells from the reference position are inside raster limits. It also display messages if some bound is out. If yo need to move the reference position you can simply write an iteration changing input position to c(i,j).

### Data input

``````# Load packages
library("raster")

# Load matrix data
m <- matrix(c(2,4,5,5,2,8,7,3,1,6,
5,7,5,7,1,6,7,2,6,3,
4,7,3,4,5,3,7,9,3,8,
9,3,6,8,3,4,7,3,7,8,
3,3,7,7,5,3,2,8,9,8,
7,6,2,6,5,2,2,7,7,7,
4,7,2,5,7,7,7,3,3,5,
7,6,7,5,9,6,5,2,3,2,
4,9,2,5,5,8,3,3,1,2,
5,2,6,5,1,5,3,7,7,2), nrow=10, ncol=10, byrow = TRUE)

# Convert matrix to RasterLayer object
r <- raster(m)

# Assign extent to raster
extent(r) <- matrix(c(0, 0, 10, 10), nrow=2)

# Plot original raster
plot(r)
text(r)
points(xFromCol(r, col=5), yFromRow(r, row=5), col="red", pch=16)
``````

### Function

``````# Function to assigning values to the adjacent raster cells based on conditions
# Input raster: RasterLayer object
# Input position: two-dimension vector (e.g. c(5,5))

AssignValuesToAdjacentRasterCells <- function(raster, position) {

# Reference position
rowPosition = position
colPosition = position

# Adjacent cells positions
adjacentBelow1 = rowPosition + 1
adjacentBelow2 = rowPosition + 2
adjacentUpper1 = rowPosition - 1
adjacentUpper2 = rowPosition - 2
adjacentLeft1 = colPosition - 1
adjacentLeft2 = colPosition - 2
adjacentRight1 = colPosition + 1
adjacentRight2 = colPosition + 2

# Check if adjacent cells positions are out of raster positions limits
belowBound1 = adjacentBelow1 <= nrow(raster)
belowBound2 = adjacentBelow2 <= nrow(raster)
upperBound1 = adjacentUpper1 > 0
upperBound2 = adjacentUpper2 > 0
leftBound1 = adjacentLeft1 > 0
leftBound2 = adjacentLeft2 > 0
rightBound1 = adjacentRight1 <= ncol(raster)
rightBound2 = adjacentRight2 <= ncol(raster)

if(upperBound2 & leftBound2) {

val1 = mean(r[adjacentUpper2:adjacentUpper1, adjacentLeft2:adjacentLeft1]) * sqrt(2)

} else {

val1 = NA

}

if(upperBound2 & leftBound1 & rightBound1) {

val2 = mean(r[adjacentUpper1:adjacentUpper2, adjacentLeft1:adjacentRight1])

} else {

val2 = NA

}

if(upperBound2 & rightBound2) {

val3 = mean(r[adjacentUpper2:adjacentUpper1, adjacentRight1:adjacentRight2]) * sqrt(2)

} else {

val3 = NA

}

if(upperBound1 & belowBound1 & leftBound2) {

val4 = mean(r[adjacentUpper1:adjacentBelow1, adjacentLeft2:adjacentLeft1])

} else {

val4 = NA

}

val5 = 0

if(upperBound1 & belowBound1 & rightBound2) {

val6 = mean(r[adjacentUpper1:adjacentBelow1, adjacentRight1:adjacentRight2])

} else {

val6 = NA

}

if(belowBound2 & leftBound2) {

val7 = mean(r[adjacentBelow1:adjacentBelow2, adjacentLeft2:adjacentLeft1]) * sqrt(2)

} else {

val7 = NA

}

if(belowBound2 & leftBound1 & rightBound1) {

val8 = mean(r[adjacentBelow1:adjacentBelow2, adjacentLeft1:adjacentRight1])

} else {

val8 = NA

}

if(belowBound2 & rightBound2) {

val9 = mean(r[adjacentBelow1:adjacentBelow2, adjacentRight1:adjacentRight2]) * sqrt(2)

} else {

val9 = NA

}

# Build matrix
mValues = matrix(data = c(val1, val2, val3,
val4, val5, val6,
val7, val8, val9), nrow = 3, ncol = 3, byrow = TRUE)

if(upperBound1) {

a = adjacentUpper1

} else {

# Warning message
cat(paste("\n Upper bound out of raster limits!"))
a = rowPosition
mValues <- mValues[-1,]

}

if(belowBound1) {

b = adjacentBelow1

} else {

# Warning message
cat(paste("\n Below bound out of raster limits!"))
b = rowPosition
mValues <- mValues[-3,]

}

if(leftBound1) {

c = adjacentLeft1

} else {

# Warning message
cat(paste("\n Left bound out of raster limits!"))
c = colPosition
mValues <- mValues[,-1]

}

if(rightBound1) {

d = adjacentRight1

} else {

# Warning
cat(paste("\n Right bound out of raster limits!"))
d = colPosition
mValues <- mValues[,-3]

}

# Convert matrix to RasterLayer object
rValues = raster(mValues)

# Assign values to raster
raster[a:b, c:d] = rValues[,]

# Assign extent to raster
extent(raster) <- matrix(c(0, 0, 10, 10), nrow = 2)

# Return raster with assigned values
return(raster)

}
``````

### Run examples

``````# Run function AssignValuesToAdjacentRasterCells

# reference position (1,1)
example1 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(1,1))

# reference position (1,5)
example2 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(1,5))

# reference position (1,10)
example3 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(1,10))

# reference position (5,1)
example4 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(5,1))

# reference position (5,5)
example5 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(5,5))

# reference position (5,10)
example6 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(5,10))

# reference position (10,1)
example7 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(10,1))

# reference position (10,5)
example8 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(10,5))

# reference position (10,10)
example9 <- AssignValuesToAdjacentRasterCells(raster = r, position = c(10,10))
``````

### Plot examples

``````# Plot examples
par(mfrow=(c(3,3)))

plot(example1, main = "Position ref. (1,1)")
text(example1)
points(xFromCol(example1, col=1), yFromRow(example1, row=1), col="red", cex=2.5, lwd=2.5)

plot(example2, main = "Position ref. (1,5)")
text(example2)
points(xFromCol(example2, col=5), yFromRow(example2, row=1), col="red", cex=2.5, lwd=2.5)

plot(example3, main = "Position ref. (1,10)")
text(example3)
points(xFromCol(example3, col=10), yFromRow(example3, row=1), col="red", cex=2.5, lwd=2.5)

plot(example4, main = "Position ref. (5,1)")
text(example4)
points(xFromCol(example4, col=1), yFromRow(example4, row=5), col="red", cex=2.5, lwd=2.5)

plot(example5, main = "Position ref. (5,5)")
text(example5)
points(xFromCol(example5, col=5), yFromRow(example5, row=5), col="red", cex=2.5, lwd=2.5)

plot(example6, main = "Position ref. (5,10)")
text(example6)
points(xFromCol(example6, col=10), yFromRow(example6, row=5), col="red", cex=2.5, lwd=2.5)

plot(example7, main = "Position ref. (10,1)")
text(example7)
points(xFromCol(example7, col=1), yFromRow(example7, row=10), col="red", cex=2.5, lwd=2.5)

plot(example8, main = "Position ref. (10,5)")
text(example8)
points(xFromCol(example8, col=5), yFromRow(example8, row=10), col="red", cex=2.5, lwd=2.5)

plot(example9, main = "Position ref. (10,10)")
text(example9)
points(xFromCol(example9, col=10), yFromRow(example9, row=10), col="red", cex=2.5, lwd=2.5)
``````

### Figure example Note: white cells mean `NA` values

For a matrix operator on a small matrix this makes sense and is tractable. However, you may want to really rethink your logic when applying a function like this to a large raster. Conceptually, this does not really track in general application. You are talking about what has traditionally been referred to as a block statistic. However, a block statistic is by nature starting at one corner of the raster and replacing blocks of values, within a specified window size, with an operator. Normally this type of operator is for aggregating data. It would be considerably more tractable if you thought in terms of using conditions to calculate a center value of a matrix. In this way you could easily use a focal function.

Just keep in mind that the raster focal function is reading in blocks of data that represent the focal values in the defined neighborhood based on the matrix passed to the w argument. The result is a vector for each neighborhood and the result of the focal operator is assigned to just the focal cell and not the entire neighborhood. Think of it as grabbing a matrix surrounding a cell value, operating on it, assigning a new value to the cell then moving to the next cell.

If you make sure that na.rm=FALSE then the vector will always represent the exact neighborhood (ie., the same length vector) and be coerced into a matrix object that can be operated on within a function. Because of this, you can simply write a function that takes the expect vector, coerces into a matrix, applies your neighborhood notation logic and then assigns a single value as the result. This function can then be passed to the raster::focal function.

Here is what would be happening at each cell based on a simple coercion and evaluation of the focal window. The "w" object would essentially be the same matrix definition that one would pass the the w argument in focal. This is what defines the size of the subset vector in each focal evaluation.

``````w=c(5,5)
x <- runif(w*w)
x <- NA
print(x)
( x <- matrix(x, nrow=w, ncol=w) )
( se <- mean(x, na.rm=TRUE) * sqrt(2) )
ifelse( as.vector(x[(length(as.vector(x)) + 1)/2]) <= se, 1, 0)
``````

Now create a function that can be applied to focal applies the above logic. In this case you could assign the se object as the value or use it as a condition in something like "ifelse" to assign a value based on an evaluation. I am adding the ifelse statement to illustrate how one would evaluate multiple conditions of the neighborhood and apply a matrix position (neighborhood notation) condition. In this dummy function the coercion of x to a matrix is completely unnecessary and there just to illustrate how it would be done. One can apply neighborhood notation conditions directly to the vector, without matrix coercion, because the position in the vector would apply to its location in the focal window and remain fixed.

``````f.rast <- function(x, dims=c(5,5)) {
x <- matrix(x, nrow=dims, ncol=dims)
se <- mean(x, na.rm=TRUE) * sqrt(2)
ifelse( as.vector(x[(length(as.vector(x)) + 1)/2]) <= se, 1, 0)
}
``````

And apply it to a raster

``````library(raster)
r <- raster(nrows=100, ncols=100)
r[] <- runif( ncell(r) )
plot(r)

( r.class <- focal(r, w = matrix(1, nrow=w, ncol=w), fun=f.rast) )
plot(r.class)
``````

You can easily update raster values by subseting raster using [row,col] notation. Just note that row and column start from upper-left corner of the raster; r[1,1] is the upper left pixel index and r[2,1] is the one underneath r[1,1]. ``````# the function to update raster cell values
focal_raster_update <- function(r, row, col) {
# copy the raster to hold the temporary values
r_copy <- r
r_copy[row,col] <- 0
#upper left
r_copy[row-1,col-1] <- mean(r[(row-2):(row-1),(col-2):(col-1)]) * sqrt(2)
#upper mid
r_copy[row-1,col] <- mean(r[(row-2):(row-1),(col-1):(col+1)])
#upper right
r_copy[row-1,col+1] <- mean(r[(row-2):(row-1),(col+1):(col+2)]) * sqrt(2)
#left
r_copy[row,col-1] <- mean(r[(row-1):(row+1),(col-2):(col-1)])
#right
r_copy[row,col+1] <- mean(r[(row-1):(row+1),(col+1):(col+2)])
#bottom left
r_copy[row+1,col-1] <- mean(r[(row+1):(row+2),(col-2):(col-1)]) * sqrt(2)
#bottom mid
r_copy[row+1,col] <- mean(r[(row+1):(row+2),(col-1):(col+1)])
#bottom right
r_copy[row+1,col+1] <- mean(r[(row+1):(row+2),(col+1):(col+2)]) * sqrt(2)
return(r_copy)
}
col <- 5
row <- 5
r <- focal_raster_update(r,row,col)

dev.set(1)
plot(r)
text(r,digits=2)
``````