Feature request: support for full joins in st_join

[HRodenhizer]

Sometimes I want to join two spatial datasets to see which geometries are in both datasets, which geometries are only in x, and which geometries are only in y. As far as I can tell, there only seems to be support for left and inner joins currently. My workaround is a somewhat annoying and inefficient task, where it is necessary to run the intersection on the inner rows twice and then remove duplicated inner rows, which is problematic when the datasets have a large number of intersecting geometries:

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

geometry = st_sfc(
  st_point(c(0, 0)),
  st_point(c(0, 1)),
  st_point(c(1, 1))
)
data_a = st_sf(a = 'a', geometry)

geometry = st_sfc(
  st_point(c(0, 0)),
  st_point(c(1, 0)),
  st_point(c(1, 1))
)
data_b = st_sf(b = 'b', geometry)
rm(geometry)

ggplot() +
  geom_sf(
    data = data_a, 
    aes(color = 'data_a'),
    alpha = 0.5
  ) +
  geom_sf(
    data = data_b, 
    aes(color = 'data_b'),
    alpha = 0.5
  )

# get all rows of a
join_a_left = data_a |>
  st_join(data_b, left = TRUE)

# get all rows of b
join_b_left = data_b |>
  st_join(data_a, left = TRUE)

# row bind join results
a_b_join = join_a_left |>
  bind_rows(join_b_left) |>
  distinct() |> # remove duplicated rows that exist in both a and b
  mutate(
    dataset = case_when(
      a == 'a' & is.na(b) ~ 'a',
      a == 'a' & b == 'b' ~ 'both',
      is.na(a) & b == 'b' ~ 'b'
    )
  )

ggplot() +
  geom_sf(
    data = a_b_join,
    aes(color = dataset)
  )

Instead, I would love to see support for full joins within st_join(). I could imagine this looking something like:

a_b_join = st_join(x = data_a, y = data_b, style = 'full')

where the argument style would replace left and could accept several options (at least "left", "inner", and "full". maybe "outer"?). I haven't really thought about the name of the argument, "style" was just the first vaguely relevant sounding term I came up with...