2

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

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()，我会收到此错误

错误： ！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")

多多边形尝试

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

4

1 回答 1

0

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

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

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

希望这可以帮助。我很高兴这些建议会让你失望。

代表

  • 您的数据
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 包于 2022-03-05 创建(v2.0.1)

于 2022-03-05T13:36:49.910 回答