m[is.na(m) & !(cbind(is.na(m[,-1L]),T) & cbind(T,is.na(m[,-ncol(m)])) & rbind(is.na(m[-1L,]),T) & rbind(T,is.na(m[-nrow(m),])))] <- 0; m; ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] ## [1,] NA NA NA NA NA NA 0 NA NA NA ## [2,] NA NA NA NA NA 0 1 0 NA NA ## [3,] NA NA NA NA NA NA 0 NA NA NA ## [4,] NA NA NA NA NA NA NA NA NA NA ## [5,] 0 NA NA NA NA NA NA NA 0 NA ## [6,] 1 0 0 0 NA 0 NA 0 1 0 ## [7,] 0 0 1 1 0 1 0 NA 0 1 ## [8,] NA NA 0 0 0 1 0 NA 0 0 ## [9,] NA 0 NA NA NA 0 NA 0 1 0 ## [10,] 0 1 0 NA NA NA NA NA 0 NA
The solution works as follows.
We construct a logical index matrix with TRUE , where the element NA And is located next to (above, below, left, or right) at least one element other than NA. Then we can index m using the logical index matrix and assign the required replacement value.
LHS logical conjunction is easy; it's just is.na(m) .
RHS logical conjunction is the hardest part. We need to perform 4 tests, one for each adjacency direction. General algorithm:
1: Index from a singular adjacency direction index that is not adjacent to any other index relative to this adjacency direction. For example, for the “right direction”, we index the leftmost column because it is not to the right of any other index. In other words, there is no column with the leftmost column on the right, so we can ignore it (and must delete it) to calculate the “right direction”.
2: Check the submatrix for NA with is.na() .
3: Then we need to bind ( cbind() to horizontal adjacent directions, rbind() for the vertical) TRUE on the opposite side (i.e. opposite the index that was deleted in step 1) of the resulting logical submatrix. This actually leads to the fact that the last index in the adjacency direction always has (pseudo) NA in its adjacent direction, so it will never be replaced due to this adjacent direction.
4: Logical AND 4 tests. The result will be a logical matrix with TRUE for elements that have NA in each adjacent cell.
5: Cancel the result of step 4. This will create a logical matrix with TRUE for elements that have at least one non-NA in any adjacent cell.
Note that there is another way to do this, which is perhaps a little intuitive. We can write each of 4 tests to test non-NA, not NA, and then logically OR them together. For the last index, you also need to bind FALSE instead of TRUE . It will look like this:
m[is.na(m) & (cbind(!is.na(m[,-1L]),F) | cbind(F,!is.na(m[,-ncol(m)])) | rbind(!is.na(m[-1L,]),F) | rbind(F,!is.na(m[-nrow(m),])))] <- 0; m;
The first approach is preferable because it requires only one negation, while the second approach requires 4 negations.
Benchmarking
library(raster); library(microbenchmark); bgoldst1 <- function(m) { m[is.na(m) & !(cbind(is.na(m[,-1L]),T) & cbind(T,is.na(m[,-ncol(m)])) & rbind(is.na(m[-1L,]),T) & rbind(T,is.na(m[-nrow(m),])))] <- 0; m; }; bgoldst2 <- function(m) { m[is.na(m) & (cbind(!is.na(m[,-1L]),F) | cbind(F,!is.na(m[,-ncol(m)])) | rbind(!is.na(m[-1L,]),F) | rbind(F,!is.na(m[-nrow(m),])))] <- 0; m; }; geotheory <- function(m) { r <- raster(m,crs='+init=epsg:27700'); extent(r) <- extent(1,ncol(m),1,nrow(m)); b <- as.matrix(buffer(r,1)); m[is.na(m) & !is.na(b)] <- 0; m; }; set.seed(1L); m <- matrix(rep(NA,100),nrow=10L); m[sample(1:100,10L)] <- 1; expected <- bgoldst1(m); identical(expected,bgoldst2(m)); ## [1] TRUE identical(expected,geotheory(m)); ## [1] TRUE microbenchmark(bgoldst1(m),bgoldst2(m),geotheory(m)); ## Unit: microseconds ## expr min lq mean median uq max neval ## bgoldst1(m) 89.380 96.0085 110.0142 107.9825 119.1015 197.149 100 ## bgoldst2(m) 87.242 97.5055 111.4725 107.3410 121.2410 176.194 100 ## geotheory(m) 5010.376 5519.7095 6017.3685 5824.4115 6289.9115 9013.201 100
set.seed(1L); NR <- 100L; NC <- 100L; probNA <- 0.9; m <- matrix(sample(c(1,NA),NR*NC,T,c(1-probNA,probNA)),NR); expected <- bgoldst1(m); identical(expected,bgoldst2(m)); ## [1] TRUE identical(expected,geotheory(m)); ## [1] TRUE microbenchmark(bgoldst1(m),bgoldst2(m),geotheory(m)); ## Unit: milliseconds ## expr min lq mean median uq max neval ## bgoldst1(m) 6.815069 7.053484 7.265562 7.100954 7.220269 8.930236 100 ## bgoldst2(m) 6.920270 7.071018 7.381712 7.127683 7.217275 16.034825 100 ## geotheory(m) 56.505277 57.989872 66.803291 58.494288 59.451588 571.142534 100