合并Shapefile和数据框

Question

我正在使用R，处理一个常规数据框(df)和一个形状文件(map2)，它们共享一个名为CD116FP的公共列。其中df有103552行，而map2只有444行。我按照以下方式加载形状文件:

map2 <- read_sf("D:/Data/tl_2019_us_cd116.shp")

我的最终目标是使用函数mapview()查看包含在map2中的地图，并显示在df中描述为"强度(intensity)"的内容，在列np_scores下。因此，我不想在map2中显示df中没有出现过的观测值。

以下是我的想法和失败经历:

1. 如果这两个对象都是常规数据框，合并它们的合理候选方法是使用merge()，但是在这种情况下应用该函数会导致结果对象失去空间属性，mapview不知道如何读取它。
2. 我尝试的另一种方法是尝试以下代码行:

map2m <- data.frame(map2, df[match(map2$CD116FP, df$CD116FP),])

但结果的维数太大(比444行要大得多)，因此当尝试绘制所需地图时，mapview会崩溃。

3. 最后，我采用了全力攻坚的方法，只是构建了一个循环来将列np添加到map2中:

map2$np=10

for (i in c(1:nrow(map2)))
{  
for (j in c(1:nrow(df)))
 {
if (identical(map2$CD116FP[i],df$CD116FP[j]))
{map2$np[i]=df$np_score[j]}
else {map2$np[i]=0}  
}
}  

然而，鉴于我的数据框的维度，这种方法需要花费太多时间。你有什么建议吗？

- Weierstraß Ramirez

---

你能分享一下 .shp 文件吗？或者分享一个类似的文件？ - Werner Hertzog

当然！在这里：https://www2.census.gov/geo/tiger/TIGER2019/CD/tl_2019_us_cd116.zip - Weierstraß Ramirez

Answer 1

我对你的数据结构有些困惑。你的df有超过100,000行，因此我猜测同一CD116FP在df中会出现多次，并且npscore可能会在这些实例之间变化。如果你想将它们合并到map2中，你需要先对它们进行聚合。

让我们尝试重建类似的设置：

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

map2 <- read_sf("C:/users/administrator/documents/shape/tl_2019_us_cd116.shp")

set.seed(69)

df <- data.frame(CD116FP = sprintf("%02d", sample(0:99, 103552, TRUE)),
                 npscores = runif(103552))

head(df)
#>   CD116FP  npscores
#> 1      95 0.6927742
#> 2      80 0.8543845
#> 3      90 0.5220353
#> 4      01 0.1449647
#> 5      76 0.9876543
#> 6      38 0.5629950

我已经让df有了与您的数据相同的行数，以展示这个解决方案可以扩展到您的问题。

让我们使用dplyr来聚合npscores：

library(dplyr)
df_sum <- df %>% 
  filter(CD116FP %in% map2$CD116FP) %>%
  group_by(CD116FP) %>%
  summarise(npscores = mean(npscores))

map2$npscores <- df_sum$npscores[match(map2$CD116FP, df_sum$CD116FP)]

现在map2已经聚合了npscores我们可以绘制图表-例如，在ggplot中：

library(ggplot2)

ggplot(map2) + 
  geom_sf(aes(fill = npscores)) +
  coord_sf(xlim = c(-180, -60),
            ylim = c(15, 70)) +
  scale_fill_gradient(low = "red", high = "gold")

在地图视图中：

library(mapview)
mapView(map2, zcol = "npscores")

^{此示例创建于2020年9月19日，使用reprex包 (v0.3.0)}

- Allan Cameron

---

谢谢！dplyr的正确使用是我所需要的。另外一个问题，当在ggplot中绘图时，如何处理npsocres中可能存在的NA值？ - Weierstraß Ramirez

实际上，这似乎不是NAs的问题。我使用mapView得到了正确的图形，但是使用ggplot时出现以下消息：“为每个多边形定义区域 Error in FUN(X[[i]], ...) : object 'np_score' not found” - Weierstraß Ramirez

我在我的代码中使用了 npscores - 只需检查名称是否一致 @WeierstraßRamirez - Allan Cameron

是的，我已经注意到了。我相应地调整了我的代码，但仍然出现了同样的错误。 - Weierstraß Ramirez

@WeierstraßRamirez，如果您的绘图代码包含名为np_score的变量，则会出现该错误。它是从哪里来的？map2中的列名是什么？ - Allan Cameron

我将你的npscores更改为np_score，以便在传递给我的时候兼容。 - Weierstraß Ramirez

Answer 2

我曾经试过在基本包中使用普通的merge，并且有了些许成功。以下是我自己工作的摘录，如果对您有价值，希望能对您有所帮助 :-)

my_data <- read_excel("TraderDataRaw.xlsx", 
      sheet = "fsa", 
      col_types= c("text","text","text","logical","numeric","numeric")) %>% 
      mutate(resp_rate=mailed/responses)

my_map <- st_read("lfsa000b16a_e.shp", stringsAsFactors = FALSE) 

my_merged_data <- merge(my_map, my_data, 
      by.x=c("CFSAUID","PRUID","PRNAME"), 
      by.y=c("CFSAUID","PRUID","PRNAME"))

- Alexa