在R中实现QGIS Dissolve功能（保留不相交要素）

我一直在使用QGIS 3.28.9-Firenze版本的Dissolve工具，该工具的*「Keep disjoint features separate（保留不相交要素）」*选项非常实用。我希望在R语言中实现该功能，无需在QGIS与R之间来回切换。

我知道R中可按分组进行融合，但输出要素数量仅与分组数一致。如下例中，6个空间要素分属box 1和box 2两组，按组融合后仅得到2个要素，未考虑要素是否共享几何关系。期望结果应为4个要素：右上角融合后的box 1、中间融合后的box 2、中间单个box 1、左下角单个box 2。

library(sf)
library(dplyr)

# Create example data
sq = function(pt, sz = 1) st_polygon(list(rbind(c(pt - sz), c(pt[1] + sz, pt[2] - sz), c(pt + sz), c(pt[1] - sz, pt[2] + sz), c(pt - sz))))
x = st_sf(box = 1:2, st_sfc(sq(c(4.2,4.2)), sq(c(0,0)), sq(c(1, -0.8)), sq(c(0.5, 1.7)), sq(c(3,3)), sq(c(-3, -3))))

# Visualise
plot(x)

# Dissolve
dissolve <- x |>
 group_by(box) |>
 summarise()
## Only two features are kept; one for box 1 and one for box 2

# Show dissolve results
dissolve
plot(dissolve)

2023年12月18日更新 #

已采纳的解决方案适用于上述多边形示例，但似乎无法支持LINESTRING（线要素）。以下是另一个示例，其中veg_type为线要素的分组列：

lines_sf <- structure(list(veg_type = c(1, 1, 4, 1, 4, 1, 1, 1), geometry = structure(list(
    structure(c(256467.281296153, 256461.531311035, 526555.12499461, 
    526505.875122078), dim = c(2L, 2L), class = c("XY", "LINESTRING", 
    "sfg")), structure(c(256461.531297717, 256442.203125007, 
    526505.875014715, 526348.312683105), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(256442.203125007, 256428.421875011, 
    526348.312500004, 526230.00012207), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(256417.468689889, 256433.437487505, 
    526217.562507417, 526347.000020735), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(256417.468689889, 256433.437487505, 
    526217.562507417, 526347.000020735), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(256433.437500004, 256457.578308105, 
    526347.000122074, 526557.687500007), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(256457.578125004, 256479.218676165, 
    526557.687500007, 526754.375008311), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(256467.281311039, 256492.01568513, 
    526555.125122067, 526754.562492598), dim = c(2L, 2L), class = c("XY", 
    "LINESTRING", "sfg"))), n_empty = 0L, crs = structure(list(
    input = "Amersfoort / RD New", wkt = "PROJCRS[\"Amersfoort / RD New\",\n    BASEGEOGCRS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4289]],\n    CONVERSION[\"RD New\",\n        METHOD[\"Oblique Stereographic\",\n            ID[\"EPSG\",9809]],\n        PARAMETER[\"Latitude of natural origin\",52.1561605555556,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",5.38763888888889,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9999079,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",155000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",463000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (X)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (Y)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone.\"],\n        BBOX[50.75,3.2,53.7,7.22]],\n    ID[\"EPSG\",28992]]"), class = "crs"), class = c("sfc_LINESTRING", 
"sfc"), precision = 0, bbox = structure(c(xmin = 256417.468689889, 
ymin = 526217.562507417, xmax = 256492.01568513, ymax = 526754.562492598
), class = "bbox"))), row.names = c(NA, -8L), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(bermtype = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"))

内容的提问来源于stack exchange，提问作者Nick