2014-04-07 119 views
1

我试图根据一个国家的省份是否与所关注的国家相邻而对空间多边形数据框进行子集划分。邻居中断

因此,假设我对捷克共和国(CZE)感兴趣,我的分析单位是第一级行政区划(请参阅情节,CZE以蓝色表示)。

如何将数据分组以包括所有捷克省以及其他国家的邻近省份,例如巴伐利亚州?

code for data frame and plot

enter image description here

回答

2

有趣的问题,这是我想出了利用poly2nb,基于this。之前我没有用过这个,所以它可能有点笨拙,但至少它似乎有效,并可能让你走上正轨。另一个选项在this线程中给出,使用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')) 

为了完整起见,在rgeos包使用gRelate替代:

# 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) 

enter image description here

+0

这很好地工作。谢谢! – BlankUsername