I'm trying to subset a spatial polygons data frame conditional on whether provinces in one country border the country of interest.
So let's say that I'm interested in the Czech Republic (CZE) and my unit of analysis are the first level administrative divisions (see plot, CZE in blue).
How could I subset the data to include all Czech provinces as well as the bordering provinces of other countries, like Bavaria for instance?
Interesting question, this is what I came up with using poly2nb, based on this. I have not worked with this before, so it might be a bit clumsy, but at least it seems to work and might get you on your way. Another option is given in this thread, using gRelate.
library(spdep)
countries<-c("CZE","DEU","AUT","SVK","POL")
EUR<-getCountries(countries,level=1) # Level 1
CZE<-subset(EUR,ISO=="CZE")
plot(EUR)
# create neighbour structure with ID_1 as names
# each item in this list contains id of all neighbouring polygons
nb = poly2nb(pl=EUR, row.names=EUR$ID_1)
# plot nb structure, collection of lines and points
plot(nb, coordinates(EUR), col=2, lty=2, add=T)
# get ID_1 names from nb items, stored in attributes
nb_attr = attributes(nb)$region.id
# take only those nb items that belong to the CZE polygons
nb1 = which(nb_attr %in% CZE$ID_1)
# convert the numeric id in the selected items to ID_1 format
neighbour_ids = unique(nb_attr[unlist(nb[nb1])])
# select polygons in EUR dataset that correspond with the found ID_1's
CZE_neighbours = subset(EUR,ID_1 %in% neighbour_ids)
# plot neighbours and CZE
plot(CZE_neighbours, col=3, add=T)
plot(CZE, col='steelblue4', add=T)
# add legend
legend('topleft',fill=c('steelblue4','green'),legend=c('CZE','neighbours'))
For completeness, an alternative using gRelate in the rgeos package:
# get relations between EUR and CZE polygons
x<-gRelate(EUR, CZE, byid = TRUE)
# Select polygon ids which border at least 1 of the CZE polygons
# ("FF2F11212" are bordering polygons)
id_bordering = EUR$ID_1[apply(x,2,function(y) sum(which(y=="FF2F11212"))>=1)]
poly_bordering = subset(EUR, ID_1 %in% id_bordering)
# plot results
plot(EUR)
plot(poly_bordering, add=T, col=4)
