R：如何计数两者之间

Question

的多边形，我试图重新显示有多少市政债券是你远离克拉科夫地图：

和城市克拉科夫改变弗罗茨瓦夫。该地图在GIMP完成。

我得到了一个shapefile文件（可在这里找到：http://www.gis-support.pl/downloads/powiaty.zip）。我读了shape文件文档包像maptools，rgdal或sf，但我无法找到一个自动功能来算，因为我不喜欢做手工。

有没有这样的功能？

业绩：地图由休伯特Szotek做了https://www.facebook.com/groups/mapawka/permalink/1850973851886654/

Answer 1

我不认为在经历了网络分析，所以我必须承认不明白的每一行代码如下所示。但它的作品！很多材料的改编从这里：https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html

这是最后的结果：

代码

# Load packages 
library(raster) # loads shapefile 
library(igraph) # build network 
library(spdep) # builds network 
library(RColorBrewer) # for plot colour palette 
library(ggplot2) # plots results 

# Load Data 
powiaty <- shapefile("powiaty/powiaty") 

首先在poly2nb函数用于计算邻近地区：

# Find neighbouring areas 
nb_q <- poly2nb(powiaty) 

这就造成我们的空间网格，我们可以在这里看到：

# Plot original results 
coords <- coordinates(powiaty) 
plot(powiaty) 
plot(nb_q, coords, col="grey", add = TRUE) 

这是我不是100％肯定发生了什么位。基本上，它是工作了所有的shape文件之间的最短距离在网络中，并返回这些对的矩阵。

# Sparse matrix 
nb_B <- nb2listw(nb_q, style="B", zero.policy=TRUE) 
B <- as(nb_B, "symmetricMatrix") 

# Calculate shortest distance 
g1 <- graph.adjacency(B, mode="undirected") 
dg1 <- diameter(g1) 
sp_mat <- shortest.paths(g1) 

已经进行的计算，现在可以将数据格式化进入绘图格式，所以最短路径矩阵合并与空间数据帧。

我不知道什么是最好的ID使用，用于参照的数据集，所以我选择了jpt_kod_je变量。

# Name used to identify data 
referenceCol <- powiaty$jpt_kod_je 

# Rename spatial matrix 
sp_mat2 <- as.data.frame(sp_mat) 
sp_mat2$id <- rownames(powiaty@data) 
names(sp_mat2) <- paste0("Ref", referenceCol) 

# Add distance to shapefile data 
powiaty@data <- cbind(powiaty@data, sp_mat2) 
powiaty@data$id <- rownames(powiaty@data) 

数据现在以适当的格式显示。使用基本功能spplot我们可以得到一个图形相当迅速：

displaylayer <- "Ref1261" # id for Krakow 

# Plot the results as a basic spplot 
spplot(powiaty, displaylayer) 

我喜欢ggplot密谋更复杂的图表，你可以控制的造型更容易。然而，它是一个比较挑剔的数据如何给进，所以我们需要重新为它的数据，我们建立图前：

# Or if you want to do it in ggplot 

filtered <- data.frame(id = sp_mat2[,ncol(sp_mat2)], dist = sp_mat2[[displaylayer]]) 
ggplot_powiaty$dist == 0 

ggplot_powiaty <- powiaty %>% fortify() 
ggplot_powiaty <- merge(x = ggplot_powiaty, y = filtered, by = "id") 
names(ggplot_powiaty) 

和剧情。我已经通过删除不需要的元素并添加了背景来定制了它。此外，为了使搜索区域位于中心黑色，我使用ggplot_powiaty[ggplot_powiaty$dist == 0, ]对数据进行子集分类，然后将其作为另一个多边形进行绘制。

ggplot(ggplot_powiaty, aes(x = long, y = lat, group = group, fill = dist)) + 
    geom_polygon(colour = "black") + 
    geom_polygon(data =ggplot_powiaty[ggplot_powiaty$dist == 0, ], 
       fill = "grey60") + 
    labs(title = "Distance of Counties from Krakow", caption = "Mikey Harper") + 
    scale_fill_gradient2(low = "#d73027", mid = "#fee08b", high = "#1a9850", midpoint = 10) + 
    theme(
    axis.line = element_blank(), 
    axis.text.x = element_blank(), 
    axis.text.y = element_blank(), 
    axis.ticks = element_blank(), 
    axis.title.x = element_blank(), 
    axis.title.y = element_blank(), 
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_blank(), 
    plot.background = element_rect(fill = "#f5f5f2", color = NA), 
    panel.background = element_rect(fill = "#f5f5f2", color = NA), 
    legend.background = element_rect(fill = "#f5f5f2", color = NA), 
    panel.border = element_blank()) 

要绘制的弗罗茨瓦夫如在文章顶部显示，只要改变displaylayer <- "Ref0264"和更新的称号。