1

I am trying to do some spatial generalised linear mixed models for areal unit data with the CARBayes package. As part of this, I need to identify the polygon neighbours into a neighbour list using spdep::poly2nb.

The shapefile I have of NSW is located here: Dropbox with shapefile

There is a freely available Australia wide version at: ABS.gov.au but I'm only after NSW.

library(spdep)
library(sf)
library(sp)
library(lwgeom)

p <- read_sf("NSW_POA_arc.shp")

p.nb <- poly2nb(p)

coords <- st_coordinates(st_centroid(st_geometry(p))) # this seems to miss some polygons

plot(st_geometry(p), border="grey")
plot(p.nb, coords, pch = 19, cex = 0.4, add=TRUE)

enter image description here

For some reason poly2nb isn't recognising that some of the polygons are neighbours to each other. I've tried using the snap feature but this only vaguely works at snap = 0.0001. It also adds connections to non-neighbours.

xx <- poly2nb(p, snap = 0.0001)
plot(st_geometry(p), border="grey")
plot(xx, coords, pch = 19, cex = 0.4, add=TRUE)

enter image description here

I'm not sure if its a topology issue but lwgeom::st_is_valid says that the polygons are valid.

Any ideas would be amazing. I'm very new to writting these questions properly so any feedback is welcome. Many thanks

  • apart from Library instead of library this is a perfectly well-written question. Clear statement, given the data and code so we can reproduce it. I subsetted it down to three of the top left regions and despite the st_distance being zero they still don't show up as connected in poly2nb... working on it... – Spacedman Mar 27 at 11:46
2

Although I thought poly2nb used the "overlaps" test for determining adjacency, I do get a different set of adjacencies if I explicitly use st_overlaps and build an nb object. This looks closer to the correct form of adjacency.

 nsw = st_read("./NSW_POA_arc.shp",quiet=TRUE)
 overlapmat = st_overlaps(nsw,sparse=FALSE)
 ## although coordinates are longitude/latitude, st_overlaps assumes that they are 
 ## planar
 ovnb = mat2listw(overlapmat)
 plot(nsw$geom);plot(ovnb, st_coordinates(st_centroid(nsw)),add=TRUE)

Produces:

enter image description here

Note this is a listw object, but the $neighbours element is a nb object.

There's some difference between this and your "extra tolerance" version:

> ovnb$neighbours
Neighbour list object:
Number of regions: 622 
Number of nonzero links: 3240 
Percentage nonzero weights: 0.8374603 
Average number of links: 5.209003 
1 region with no links:
132
> nbsnap <- poly2nb(nsw, snap = 0.0001)
> nbsnap
Neighbour list object:
Number of regions: 622 
Number of nonzero links: 3260 
Percentage nonzero weights: 0.8426298 
Average number of links: 5.241158 
1 region with no links:
132

For example the second region has two more neighbours in the extra tolerance version than in the overlaps version:

> ovnb$neighbours[[2]]
[1]  99 155 159 160
> nbsnap[[2]]
[1]  74  94  99 155 159 160

enter image description here

Here region 2 is in yellow, the four greenish ones are the "overlaps" regions, the two blueish ones are the additional ones found by increasing the tolerance.

I don't see any evidence of the increased snap tolerance adding neighbours where there's no close connection, except as above where it has added a connection where there seems to be a single point common to the yellow and one of the blueish polygons.

Another approach to adjacency is to compute st_distance and apply a near-zero cutoff. This should work well and is flexible enough that you can increase the tolerance enough to jump over digitising errors or very small gaps caused by narrow rivers between regions.

| improve this answer | |
  • Hi @Spacedman, thank you for your answer. Do you know of a way to individually assign neighbours to a neighbour list? If that is possible to upload a csv or similar with them pre-populated? If both of those are possible is there a way to have the map interactive with labels (leaflet?) and then populate the list from there. Or is all of that a bit much because it won't affect the CARBayes calculations. – hclutter Mar 30 at 5:12
  • Yes, there's a whole bunch of functions in spdep for converting, creating, and modifying listw and nb objects, including an interactive (but clunky) edit function. Underneath it all an nb is a list and you can mostly treat it as such by setting elements with mynb[[3]]=c(4,5,7) – Spacedman Mar 30 at 6:48
  • Thanks for your help. I'll give that a go and see how I go. – hclutter Mar 31 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.