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

Question

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

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 session 中所有打开的文档并关闭 R。重新打开新的 R session 并仅打开包含要测试的代码的文件。
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 package 于 2022 年 3 月 5 日创建(v2.0.1)}