R gis: add polygon attribute based on neighbors value?

Question

I have a set of forest stands (SpatialPolygonDataFrame) distributed randomly over the landscape, i.e. dispersed and aggregated. For each polygon, I want to decide if it has an open edge or not. Polygon has open edge if:

• has no neighbors
• has no neighbors at least on one side;
• has neighbors, but the difference between the tree height between stand and neighbors is more than 5

I wander how to add attribute open_edge = TRUE/FALSE to individual polygons? In raster package, there is a promising approach using moving window. However, my original data are feature classes, and unfortunately are not so raster like as in the working example.

I though to (pseudocode):

• subset each stand one by one (in for loop)
• create surrounding buffer
• subset surrounding stand by overlap of buffer with the stands
• if any neighbors -> compare the height. If difference > 5, open_edge = TRUE

But, this approach does not consider what the stand has let's say neighbors only at 3 sides, i.e. as rook neighborhood. The poly2nb tool seems promising, but how to add attributes to individual stands?

Here is my dummy approach but I wonder if you have a more efficient solution?

Create dummy data:

library(ggplot2)  # for choropleth map plot
library(broom) # to convert spatial data to dataframe
library(mapproj)
library(spdep)    # neighbours
library(rgdal)
library(rgeos)
library(sf)
library(raster)
library(dplyr)
library(spData)
library(sf)

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, 1, 1, 1,
                             NA, NA, NA, 1, 9, 1,
                             NA, NA, NA, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Convert raster to polygon
polys <- rasterToPolygons(r)

Identify if stand has open edge, on example of one stand:

# Subset first row in SpatialPolygonDataFrame
i = 10
one = polys[i, ]

# Keep the remaining polygons
left = polys[-i,]

# Create buffer within distance
buff = buffer(one, width = 100)

# subset set of neighbours by spatial overlap
nbrs <- left[which(gContains(sp::geometry(buff),
                                    sp::geometry(left), byid = T)),    
# Compare if the values are different
height.one  = rep(one$layer[1], nrow(nbrs))
height.nbrs = nbrs$layer

# Get the differences between the neighbouring stands
difference = height.one - height.nbrs

# If the difference in at least one stand is 
# in more than 5, set open_edge = TRUE 
# or if no neighbours find
one$open_edge <- any(difference > 5)

Answer 1

It seems like ther is a more simple solution. You can use focal a sliding window function from raster package.

Here is an example:

library(raster)

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, 2, 1, 3, 1, 
                             NA, NA, 1, 1, 1, 1,
                             NA, NA, 1, 2, 2, 1,
                             NA, NA, 1, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Prepare function for sliding window
focal_func <- function(x) {
    # Assuming 3x3 moving window
    # central cell has index 5
    # Check if the cell is not NA
    if (is.na(x[5])){
        return(NA)
    }

    # Check if any surrounding is NA
    if (any(is.na(x[-5]))) {
        return(TRUE)
    }

    # Check difference
    if (any((x[5] - x[-5]) > 5)) {
        return(TRUE)
    }

    # Else, return FALSE
    return(FALSE)
}

# Apply focal_function with sliding window
res <- focal(r, w = matrix(rep(1, 9), 3), fun = focal_func, pad = TRUE)

# Check if the same as desired output
res_mat <- as.matrix(res)
res[!is.na(res)] == 1

It gives:

[1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE
[13]  TRUE  TRUE  TRUE  TRUE  TRUE

I.e. the same as desired output.

Answer 2

Based on @RobertHijmans answer and how to get list of neighbors I have created a step by step approach to check my neighbors and evaluate height difference.

Step by step:

1. checks how many neighbors the cell has If number of neighbors is < 8, set the open_edge as TRUE, otherwise FALSE
2. Loop through index of cells that have 8 neighbors (queen case): compare the height difference between central and surrounding cells. If any of the difference is > 5, than open_edge is set to TRUE.

Here is slightly more complicated situation, allowing more stands to have a neighbors:

Create dummy data:

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, 2, 1, 3, 1, 
                             NA, NA, 1, 1, 1, 1,
                             NA, NA, 1, 2, 2, 1,
                             NA, NA, 1, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Convert raster to polygon
polys <- rasterToPolygons(r)

Calculate neighbors based on continuity, taking into account possible gaps or slivers between stands:

# Calculate the count and position of neighbors; this allows to create one by one list
nb <- poly2nb(polys, 
              snap = 10) # snap corrects for the gaps/slivers

# store the number of neighbors for each central cell
polys$nb_count<- card(nb)

# Has the stand an open edge? Is surrounded by neighbors, pre-value is FALSE
polys$open_edge = ifelse(card(nb) <8, "TRUE", "FALSE")

If a cell has complete neighbors, compare the differences in height between them:

# Get the position of the cell surrounded by neighbors
center.index <- which(polys$nb_count == 8)

# Get the stand height of a stand
# as a vector to compare element wise
center.height = polys[center.index,]$layer

# Loop through the cells with neighbors:
# - keep height of the central stand
# - get height of neighbors
# compare the height between them
# if difference is more than 5: => open_edge = TRUE
for (i in seq_along(center.index)) {

  # Get central stand height 
  center.height = polys[center.index[i],]$layer

  # Identify neighbors of the stands
  # by the index value
  nb.index = unlist(nb[center.index[i]])

  # Get heights of the stands
  nb.height = polys[nb.index,]$layer

  # Adjust Center.height length as a vector to allow element wise comparison
  center.height.v = rep(center.height, length(nb.index))

  # Compare the heights 
  h.diff = center.height.v - nb.height

  # if any difference is more than 5 => change open_edge = TRUE
  if (any(h.diff > 5)) {
    polys@data[center.index[i],]$open_edge <- "TRUE"
  }
}

Look at final data output:

> polys@data
   layer nb_count open_edge
1      9        0      TRUE
2      2        3      TRUE
3      1        5      TRUE
4      3        5      TRUE
5      1        3      TRUE
6      1        5      TRUE
7      1        8     FALSE
8      1        8     FALSE
9      1        5      TRUE
10     1        5      TRUE
11     2        8     FALSE
12     2        8     FALSE
13     1        5      TRUE
14     1        3      TRUE
15     1        5      TRUE
16     1        5      TRUE
17     1        3      TRUE

Answer 3

To get you started with spdep::poly2nb

library(raster)

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, 1, 1, 1,
                             NA, NA, NA, 1, 9, 1,
                             NA, NA, NA, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Convert raster to polygon
polys <- rasterToPolygons(r)

library(spdep)
nb <- poly2nb(polys)

plot(polys)
text(polys, 1:length(polys))

str(nb)
#List of 10
# $ : int 0
# $ : int [1:3] 3 5 6
# $ : int [1:5] 2 4 5 6 7
# $ : int [1:3] 3 6 7
# ...

So poly 1 has no neighbors, poly 2 has neighbors 3, 5, 6, etc.

Now you can use sapply to compute things. For example the number of neighbors

nbcnt <- sapply(nb, function(i) length(i[i>0]))
nbcnt 
#[1] 0 3 5 3 5 8 5 3 5 3

To add this back to the polygons

polys$nbcnt <- nbcnt