dplyr full_join on geometry columns of sf objects

Question

I am trying to perform a full join with two sf objects, based on equality of their geometries. I know the two objects have the same geometry type and the same crs.

The st_join function for spatial joins can use the st_equals predicate but only allows left joins and inner joins (for good reasons, I assume). However, when converting the sf objects to tibbles (while keeping the geometry list column present), I can perform a regular full_join using dplyr, and it works exactly the way I want, at least for a small toy example:

suppressPackageStartupMessages({
  library(sf)
  library(dplyr)
  library(tibble)
})

# Create some points
p1 = st_sfc(st_point(c(0,0)), crs = 4326)
p2 = st_sfc(st_point(c(1,0)), crs = 4326)
p3 = st_sfc(st_point(c(0,1)), crs = 4326)

# Create two sf objects.
# Both contain two points, of which one is shared between the two objects.
sf1 = st_as_sf(data.frame(foo = c(1, 2), geometry = c(p1, p2)))
sf2 = st_as_sf(data.frame(bar = c("a", "b"), geometry = c(p1, p3)))

sf1
#> Simple feature collection with 2 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 0 ymin: 0 xmax: 1 ymax: 0
#> CRS:            EPSG:4326
#>   foo    geometry
#> 1   1 POINT (0 0)
#> 2   2 POINT (1 0)
sf2
#> Simple feature collection with 2 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 0 ymin: 0 xmax: 0 ymax: 1
#> CRS:            EPSG:4326
#>   bar    geometry
#> 1   a POINT (0 0)
#> 2   b POINT (0 1)

# Convert to tibble for a dplyr full join on the geometry column.
tbl_joined = full_join(as_tibble(sf1), as_tibble(sf2), by = "geometry")

# Convert back to sf.
sf_joined = st_as_sf(tbl_joined)

sf_joined
#> Simple feature collection with 3 features and 2 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 0 ymin: 0 xmax: 1 ymax: 0
#> CRS:            EPSG:4326
#> # A tibble: 3 x 3
#>     foo    geometry bar  
#>     
#> 1     1       (0 0) a    
#> 2     2       (1 0)  
#> 3    NA       (0 1) b

^{Created on 2020-10-15 by the reprex package (v0.3.0)}

That is great, but it does not really feel like a "correct" way of doing this. My question therefore is: for this toy example it works fine, but could there be disadvantages of this approach (joining geometry columns without a true spatial join function) that I did not think of?

dplyr full_join on geometry columns of sf objects

0 Answers