R 凸包方法产生多个多边形

Question

我可以使用以下代码在 R 中围绕我的点制作凸包：

library(sf)
library(tidyverse)

chull_polytd <- td %>%
  st_transform(., crs="EPSG:4326") %>%
  summarise( geometry = st_combine( geometry ) ) %>%
  st_buffer(., dist = 20000) %>%
  st_concave_hull(.,
                  dfMaxLength = 50000,
                  allow_holes = FALSE,
                  ratio = 0.1)

ggplot() +
  geom_sf(data = chull_polytd)+
  geom_sf(data =td)

但我希望能够将其制作成 3 个多边形。我希望 st_convex_hull 中有一个我可以控制的参数，但我认为没有。我想知道是否有人知道如何在不首先定义组的情况下生成 3 个多边形？

以下是数据：

> dput(td)
structure(list(SNES_name = c("Heleioporus australiacus", "Heleioporus australiacus", 
"Heleioporus australiacus", "Heleioporus australiacus", "Heleioporus australiacus", 
"Heleioporus australiacus", "Heleioporus australiacus", "Heleioporus australiacus", 
"Heleioporus australiacus", "Heleioporus australiacus", "Heleioporus australiacus", 
"Heleioporus australiacus", "Heleioporus australiacus"), geometry = structure(list(
    structure(c(1242924.08515837, -3928778.41765649), class = c("XY", 
    "POINT", "sfg")), structure(c(1389765.05347474, -3919507.94301041
    ), class = c("XY", "POINT", "sfg")), structure(c(1735121.88594707, 
    -3851559.44018895), class = c("XY", "POINT", "sfg")), structure(c(1727184.40284067, 
    -3850191.62128326), class = c("XY", "POINT", "sfg")), structure(c(1727132.24919221, 
    -3850076.78690896), class = c("XY", "POINT", "sfg")), structure(c(1745144.80155694, 
    -3852642.07181568), class = c("XY", "POINT", "sfg")), structure(c(1719953.84894257, 
    -3848668.53251429), class = c("XY", "POINT", "sfg")), structure(c(1747196.75400538, 
    -3852739.36884773), class = c("XY", "POINT", "sfg")), structure(c(1730797.99218133, 
    -3849885.60286591), class = c("XY", "POINT", "sfg")), structure(c(1747249.88557053, 
    -3852321.90609061), class = c("XY", "POINT", "sfg")), structure(c(1746415.57371839, 
    -3851800.11378268), class = c("XY", "POINT", "sfg")), structure(c(1739547.5569871, 
    -3850660.90862698), class = c("XY", "POINT", "sfg")), structure(c(1747894.13898112, 
    -3851838.70406853), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 1242924.08515837, 
ymin = -3928778.41765649, xmax = 1747894.13898112, ymax = -3848668.53251429
), class = "bbox"), crs = structure(list(input = "GDA94 / Australian Albers", 
    wkt = "PROJCRS[\"GDA94 / Australian Albers\",\n    BASEGEOGCRS[\"GDA94\",\n        DATUM[\"Geocentric Datum of Australia 1994\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4283]],\n    CONVERSION[\"Australian Albers\",\n        METHOD[\"Albers Equal Area\",\n            ID[\"EPSG\",9822]],\n        PARAMETER[\"Latitude of false origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",132,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",-18,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",-36,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Statistical analysis.\"],\n        AREA[\"Australia - Australian Capital Territory; New South Wales; Northern Territory; Queensland; South Australia; Tasmania; Western Australia; Victoria.\"],\n        BBOX[-43.7,112.85,-9.86,153.69]],\n    ID[\"EPSG\",3577]]"), class = "crs"), n_empty = 0L)), row.names = c("...107", 
"...137", "...400", "...401", "...402", "...403", "...404", "...405", 
"...406", "...407", "...408", "...409", "...410"), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(SNES_name = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"))

Answer 1

您可以使用

st_cast()

:

取消汇总缓冲区的分组

library(sf)
library(dplyr)
library(ggplot2)

chull_polytd <- td %>%
  st_transform(., crs="EPSG:4326") %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_buffer(., dist = 20000) %>%
  st_cast("POLYGON") %>%
  mutate(buffid = 1:n()) # For illustrative/plot purposes only

chull_polytd
# Simple feature collection with 3 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 145.5026 ymin: -35.54839 xmax: 151.3269 ymax: -33.90024
# Geodetic CRS:  WGS 84
#                         geometry buffid
# 1 POLYGON ((150.826 -34.28015...      1
# 2 POLYGON ((147.474 -35.26087...      2
# 3 POLYGON ((145.947 -35.41275...      3

结果如下：

ggplot() +
  geom_sf(data = chull_polytd) +
  geom_sf(data = td) +
  geom_sf_text(data = chull_polytd,
               aes(label = buffid),
               size = 4,
               hjust = 6,
               fun.geometry = st_centroid,
               colour = "#920000")