在 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")

但是，如果我按照其他 sites/questions 的建议取消对 group_by() 的注释，则会收到此错误

错误：！ tibble 中的所有列都必须是向量。 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

作为对您的评论的 follow-up，我准备了一个 reprex，以便您可以测试代码。它应该工作...

如果还是不行，这里有一些建议：

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

希望这对您有所帮助。我祈祷这些建议会让你摆脱困境。

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 (v2.0.1)}

于 2022 年 3 月 5 日创建

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

Create Polygons by Category in sf file R

r

polygon

tidyverse

sf

mapview