I have been using the dissolve function in QGIS (I'm using version 3.28.9-Firenze). What I appreciate about this function is that it has an option to 'Keep disjoint features separate'. I'm looking for a way to do this in R, so that I don't need to switch to QGIS back and forward.
How would I do this in R? The example shows polygons, but I would like to do this for line features as well. Below is an example dataset and what I've tried in R.
I know dissolve by groups can be done in R, but the output only keeps the samen number of features as there are groups in the dataset. In the example below, I have 6 different spatial features, either beloning to box 1 or box 2. When I dissolve by group, the result is just two features, one for each groups. It does not take into account whether the dissolved features share geometry. The desired result should be four different features (topright - dissolved box 1; middle - dissolved box 2; middle - single box 1; bottomleft - single 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)
EDIT 18-Dec-2023: The accepted solution works for the example I provided above (polygons), but does not seem to work for LINESTRINGS. Here's a different example where 'veg_type' is the column to which group each line segment belows to. :
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"))
perhaps this wil work?
x %>% group_by(box) %>% summarise() %>% sf::st_cast("POLYGON")
# box
# <int> <POLYGON>
# 1 1 ((0 0.2, 2 0.2, 2 -1.8, 0 -1.8, 0 0.2))
# 2 1 ((2 2, 2 4, 3.2 4, 3.2 5.2, 5.2 5.2, 5.2 3.2, 4 3.2, 4 2, 2 2))
# 3 2 ((-4 -2, -2 -2, -2 -4, -4 -4, -4 -2))
# 4 2 ((-1 -1, -1 1, -0.5 1, -0.5 2.7, 1.5 2.7, 1.5 0.7, 1 0.7, 1 -1, -1 -1))