在 sf 文件 R 中按类别创建多边形

Question

我有来自关键和濒危习惯联邦登记处的大量坐标。我正在尝试将这些地图数字化以进行分析。以下是数据样本作为示例。

library(tidyverse)
library(mapview)
library(sf)
library(ggplot2)

lat <- c(300786, 301348,  301467, 302384, 304644, 304747, 304863, 304882)
lon <- c(4215918, 4215650, 4215784, 4216077, 4211303, 4211336, 4211395, 4211457)
geodeticDA <- c("NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83")
utmZone <- c("11N", "11N", "11N", "11N", "11N", "11N", "11N", "11N")
ID <- c("A", "A", "A", "A", "B", "B", "B", "B")

df <- data.frame(lat, lon, geodeticDA, utmZone, ID)

我已经成功地将 df 转换为 sf 并使用以下

plot_local1 <- st_as_sf(df, 
                        coords = c("lat", "lon"),
                         crs = "+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")

mapview(plot_local1, map.types = "Esri.NatGeoWorldMap") 

将其

绘制成图形我的主要目标是根据 ID（A 或 B）创建两个不同的多边形，但我很挣扎。

我可以用以下命令将它们全部变成一个多边形：

polygon <- plot_local1 %>% 
# dplyr::group_by(ID) %>% 
  dplyr::summarise() %>%
  st_cast("POLYGON")

mapview(polygon, map.types = "Esri.NatGeoWorldMap") 

但是，如果我按照建议取消注释 group_by()通过其他网站/问题，我收到此错误

错误： ！ tible 中的所有列都必须是向量。 x 列geometry 是一个sfc_POINT/sfc 对象。

我已经设法使用以下代码创建多个多边形，但是当我绘制它时，它不正确。

polys <- st_sf(
  aggregate(
    plot_local1$geometry,
    list(plot_local1$ID),
    function(g){
       st_cast(st_combine(g),"POLYGON")
    }
   ))

mapview(polys, map.types = "Esri.NatGeoWorldMap") 

当然，我可以为每个多边形创建单独的数据框，然后稍后将它们组合在一起，但考虑到我的数据量我正在处理这确实不切实际。我还更新了我正在使用的所有软件包，所以我知道这不是问题。感谢任何和所有的帮助！

Answer 1

作为您评论的后续行动，我准备了一个表示，以便您可以测试代码。它应该可以工作...

如果它不起作用，这里有一些建议：

1. 确保所有库都是最新的，尤其是 tidyverse、mapview 和 SF。在我这边，我使用以下版本运行代码：tidyverse 1.3.1、mapview 2.10.0 和 sf 1.0.6
2. 关闭所有打开 R 会话中的文档并关闭 R。重新打开新的 R 会话并仅打开包含要测试的代码的文件。
3. 仅加载运行代码所需的库。

希望这有帮助。我祈祷这些建议能让你摆脱困境。

Reprex

• 您的数据

library(tidyverse)
library(mapview)
library(sf)

lat <- c(300786, 301348,  301467, 302384, 304644, 304747, 304863, 304882)
lon <- c(4215918, 4215650, 4215784, 4216077, 4211303, 4211336, 4211395, 4211457)
geodeticDA <- c("NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83")
utmZone <- c("11N", "11N", "11N", "11N", "11N", "11N", "11N", "11N")
ID <- c("A", "A", "A", "A", "B", "B", "B", "B")

df <- data.frame(lat, lon, geodeticDA, utmZone, ID)


plot_local1 <- st_as_sf(df, 
                        coords = c("lat", "lon"),
                        crs = "+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")

mapview(plot_local1, map.types = "Esri.NatGeoWorldMap") 

• 应该有效的代码...

polygon <- plot_local1 %>% 
  dplyr::group_by(ID) %>% 
  dplyr::summarize() %>% 
  sf::st_cast("POLYGON")

mapview(polygon, map.types = "Esri.NatGeoWorldMap") 

^{由 reprex 包 (v2.0.1)}