Unambiguous Neighborhoods with spdep

Interactive Version

(0) Initialization Code

library(tidyverse) |> suppressPackageStartupMessages()
library(mapview) |> suppressPackageStartupMessages()
# Fixes "API Key Required" message in mapview
mapviewOptions(basemaps = c("Esri.WorldGrayCanvas","Stadia.AlidadeSmoothDark","OpenStreetMap", "Esri.WorldImagery", "OpenTopoMap"))
library(rnaturalearth) |> suppressPackageStartupMessages()
library(spdep) |> suppressPackageStartupMessages()

• Relevant Moraga
• Relevant Pebesma

(1) Finding Neighbors in Asia

(1.1) Download Countries as POLYGONs from rnaturalearth

Here we use rnaturalearth::ne_countres(), with the continent option. We specify scale=50 since otherwise smaller countries like Bahrain get excluded!

asia_countries_sf <- rnaturalearth::ne_countries(
  continent = "Asia", scale=50
)
mapview(asia_countries_sf, label="geounit")

(1.2) Detecting Neighbors with spdep::poly2nb()

Given an sf object, spdep::poly2nb() identifies all pairs of touching POLYGONs.

library(spdep)
asia_countries_nb <- asia_countries_sf |>
  spdep::poly2nb(
    queen = TRUE,
    row.names=asia_countries_sf$geounit
  )

Warning in spdep::poly2nb(asia_countries_sf, queen = TRUE, row.names = asia_countries_sf$geounit): some observations have no neighbours;
if this seems unexpected, try increasing the snap argument.

Warning in spdep::poly2nb(asia_countries_sf, queen = TRUE, row.names = asia_countries_sf$geounit): neighbour object has 9 sub-graphs;
if this sub-graph count seems unexpected, try increasing the snap argument.

asia_countries_nb

Neighbour list object:
Number of regions: 53 
Number of nonzero links: 164 
Percentage nonzero weights: 5.838377 
Average number of links: 3.09434 
7 regions with no links:
Taiwan, Sri Lanka, Singapore, Philippines, Japan, Bahrain, Indian Ocean
Territories
9 disjoint connected subgraphs

nb_cents <- sf::st_point_on_surface(asia_countries_sf) |>
  sf::st_as_sf()

Warning: st_point_on_surface assumes attributes are constant over geometries

Warning in st_point_on_surface.sfc(st_geometry(x)): st_point_on_surface may not
give correct results for longitude/latitude data

mapview(nb_cents, label="geounit")

nb_lines <- spdep::nb2lines(
  asia_countries_nb,
  coords=st_coordinates(nb_cents),
  as_sf=TRUE
)
nb_lines$edgelabel <- paste0(nb_lines$i_ID, "-", nb_lines$j_ID)
mapview(nb_cents, label="geounit") + mapview(nb_lines, label="edgelabel")

print(asia_countries_nb)

Neighbour list object:
Number of regions: 53 
Number of nonzero links: 164 
Percentage nonzero weights: 5.838377 
Average number of links: 3.09434 
7 regions with no links:
Taiwan, Sri Lanka, Singapore, Philippines, Japan, Bahrain, Indian Ocean
Territories
9 disjoint connected subgraphs

asia_nb_names <- attributes(asia_countries_nb)$region.id
afg_index <- which(asia_nb_names == "Afghanistan")
asia_countries_sf$neighbors <- "other"
asia_countries_sf$neighbors[afg_index] <- "area"
asia_countries_sf$neighbors[asia_countries_nb[[afg_index]]] <- "neighbors"
afg_sf <- asia_countries_sf |> filter(neighbors != "other")
ggplot(afg_sf) +
  geom_sf(aes(fill = neighbors)) +
  theme_bw() +
  scale_fill_manual(
    values = c("gray30", "gray", "white")
  )

mapview(afg_sf, zcol="neighbors", label="geounit")

nb_weights <- asia_countries_nb |>
  spdep::nb2listw(zero.policy=TRUE)
nb_weights

Characteristics of weights list object:
Neighbour list object:
Number of regions: 53 
Number of nonzero links: 164 
Percentage nonzero weights: 5.838377 
Average number of links: 3.09434 
7 regions with no links:
Taiwan, Sri Lanka, Singapore, Philippines, Japan, Bahrain, Indian Ocean
Territories
9 disjoint connected subgraphs

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 46 2116 46 33.52657 234.6487

afg_df <- tibble(
  nb_name=asia_countries_sf$geounit[nb_weights$neighbours[[afg_index]]],
  nb_index=nb_weights$neighbours[[afg_index]],
  nb_weight=nb_weights$weights[[afg_index]]
)
afg_df

nb_name	nb_index	nb_weight
Uzbekistan	3	0.1666667
Turkmenistan	5	0.1666667
Tajikistan	9	0.1666667
Pakistan	18	0.1666667
Iran	34	0.1666667
China	40	0.1666667

nb_dist_weights <- asia_countries_nb |>
  spdep::nb2listwdist(nb_cents, zero.policy=TRUE)
nb_dist_weights

Characteristics of weights list object:
Neighbour list object:
Number of regions: 53 
Number of nonzero links: 164 
Percentage nonzero weights: 5.838377 
Average number of links: 3.09434 
7 regions with no links:
Taiwan, Sri Lanka, Singapore, Philippines, Japan, Bahrain, Indian Ocean
Territories
9 disjoint connected subgraphs

Weights style: raw 
Weights constants summary:
     n   nn      S0       S1       S2
raw 46 2116 28.5084 27.65386 123.5441

afg_dist_df <- tibble(
  nb_name=asia_countries_sf$geounit[nb_dist_weights$neighbours[[afg_index]]],
  nb_index=nb_dist_weights$neighbours[[afg_index]],
  nb_weight=nb_dist_weights$weights[[afg_index]]
)
afg_dist_df

nb_name	nb_index	nb_weight
Uzbekistan	3	0.1306825
Turkmenistan	5	0.1251905
Tajikistan	9	0.1321491
Pakistan	18	0.1657872
Iran	34	0.0890133
China	40	0.0298321