Enable sf::st_as_sf to convert neighbors list object with class nb to sf

[alexyshr]

Dear Dr. Edzer,

Maybe it is a good idea to enable sf::st_as_sf to convert a neighbors list object with class nb (output of spdep::poly2nb) to simple features. Here I am providing two functions ( nb_to_df and st_segment ) to solve the conversion (I got those from the Internet).

Sorry, but I need to take apart some time to learn it and do it by myself.

Best,

library(spdep)
library(ggplot2)
library(sf)

# function converts nb object to a data.frame
# with the coordinates of the lines "x"    "y"    "xend" "yend"
  nb_to_df <- function(nb, coords){
    x <- coords[, 1]
    y <- coords[, 2]
    n <- length(nb)
    cardnb <- card(nb)
    i <- rep(1:n, cardnb)
    j <- unlist(nb)
    return(data.frame(x=x[i], y=y[i], 
                      xend=x[j], yend=y[j]))
  }

# Use this function to convert the dataframe to sf
st_segment = function(r){sf::st_linestring(t(matrix(unlist(r), 2, 2)))} 
    #as we are not using here byrow=T, we have to use
    #transpose to have the format:
    #x1 y1
    #X2 Y2

#Here one example
ColData1 <- sf::st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
coords1 = sf::st_coordinates(sf::st_centroid(ColData1)) # coordinates of the census tract centroids
#> Warning in st_centroid.sf(ColData1): st_centroid assumes attributes are
#> constant over geometries of x
IDs1 = row.names(ColData1)
ColData1$IDs1 = IDs1
col_nb1 = poly2nb(ColData1, queen = TRUE)
summary(col_nb1)
#> Neighbour list object:
#> Number of regions: 49 
#> Number of nonzero links: 236 
#> Percentage nonzero weights: 9.829238 
#> Average number of links: 4.816327 
#> Link number distribution:
#> 
#>  2  3  4  5  6  7  8  9 10 
#>  5  9 12  5  9  3  4  1  1 
#> 5 least connected regions:
#> 1 6 42 46 47 with 2 links
#> 1 most connected region:
#> 20 with 10 links
class(col_nb1) # check the class of col_nb
#> [1] "nb"
#A neighbours list with class nb (output of spdep::poly2nb )

# this is a dataframe with the coordinates of the lines "x"    "y"    "xend" "yend"
col_nb1_df = nb_to_df(col_nb1, coords1)

col_nb1_df$geom = st_sfc(sapply(1:nrow(col_nb1_df), 
                                function(i){st_segment(col_nb1_df[i,])}, simplify=FALSE)) 
          # the simplify = FALSE force to return LINESTRING (x1 y1, x2 y2)

#convert to sf
col_nb1_sf = sf::st_sf(col_nb1_df)

# plot
ggplot() +
  geom_sf(data=ColData1, fill=NA, color="red", lwd=0.5) +
  geom_sf(data=sf::st_centroid(ColData1), pch=19, color="blue") +
  geom_sf(data=col_nb1_sf, color="blue", lwd=0.3, lty=2) +
  geom_sf_label(data=ColData1, aes(label=IDs1), size = 2) +
  #geom_sf_text(data=ColData1, aes(label=IDs1), size = 3) +
  labs(title="Adjacency Connection", subtitle = "Columbus Shapefile") + 
  xlab("") +
  ylab("")
#> Warning in st_centroid.sf(ColData1): st_centroid assumes attributes are
#> constant over geometries of x

^{Created on 2023-03-14 with reprex v2.0.2}

[rsbivand]

Where on the internet and when? Did you miss spdep::nb2lines(nb, wts, coords, proj4string=NULL, as_sf=FALSE)? Yes, as_sf= could be TRUE. See https://github.com/r-spatial/spdep/blob/82514fb81da0c6c5285bf4219d7979dac34acec4/R/nb2lines.R#L4

[alexyshr]

I modified the function nb_to_df from here. Jun 15, 2015. I got the function st_segment from here. Feb 14, 2019.

Below is the right solution with spdep::nb2lines.

Thank you for your answer.

library(spdep)
library(ggplot2)
library(sf)

ColData1 <- sf::st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
coords1 = sf::st_coordinates(sf::st_centroid(ColData1)) # coordinates of the census tract centroids
#> Warning in st_centroid.sf(ColData1): st_centroid assumes attributes are
#> constant over geometries of x
IDs1 = row.names(ColData1)
ColData1$IDs1 = IDs1
col_nb1 = poly2nb(ColData1, queen = TRUE)
summary(col_nb1)
#> Neighbour list object:
#> Number of regions: 49 
#> Number of nonzero links: 236 
#> Percentage nonzero weights: 9.829238 
#> Average number of links: 4.816327 
#> Link number distribution:
#> 
#>  2  3  4  5  6  7  8  9 10 
#>  5  9 12  5  9  3  4  1  1 
#> 5 least connected regions:
#> 1 6 42 46 47 with 2 links
#> 1 most connected region:
#> 20 with 10 links
class(col_nb1) # check the class of col_nb
#> [1] "nb"

# convert nb to sf with binary weights
col_nb1_sf = spdep::nb2lines(col_nb1, coords=coords1, proj4string=NULL, as_sf=T) # weights are binary

# plot
ggplot() +
  geom_sf(data=ColData1, fill=NA, color="red", lwd=0.5) +
  geom_sf(data=sf::st_centroid(ColData1), pch=19, color="blue") +
  geom_sf(data=col_nb1_sf, color="blue", lwd=0.3, lty=2) +
  geom_sf_label(data=ColData1, aes(label=IDs1), size = 2) +
  #geom_sf_text(data=ColData1, aes(label=IDs1), size = 3) +
  labs(title="Adjacency Connection", subtitle = "Columbus Shapefile") + 
  xlab("") +
  ylab("")
#> Warning in st_centroid.sf(ColData1): st_centroid assumes attributes are
#> constant over geometries of x

^{Created on 2023-03-15 with reprex v2.0.2}

[edzer]

Thanks, @alexyshr !