FR: Implement `by_feature` to `st_union()` when x and y are provided

[olivroy]

Hi,

I was wondering if it would be possible to implement a method by_feature argument to st_union() when x and y a provided.

My use case would be to have two sf objects as data frame. and to compute pairwise union

I have 2 shapefiles that have the same number of rows. I want to create a line for each row.

I can't figure out how to do it.

Here is my example code

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
# For the sake of reproducibility
tab1 <- data.frame(
  long = c(-129, -80, -100),
  lat = c(56, 69, 60)
)
tab2 <- data.frame(
  long = c(-100, -90, -120),
  lat = c(59, 75, 62)
)
my_x <- sf::st_as_sf(tab1, coords = c("long", "lat"), crs = 4326)
my_y <- sf::st_as_sf(tab2, coords = c("long", "lat"), crs = 4326)

# I have to start with sf objects

my_x
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -129 ymin: 56 xmax: -80 ymax: 69
#> Geodetic CRS:  WGS 84
#>          geometry
#> 1 POINT (-129 56)
#> 2  POINT (-80 69)
#> 3 POINT (-100 60)
my_y
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -120 ymin: 59 xmax: -90 ymax: 75
#> Geodetic CRS:  WGS 84
#>          geometry
#> 1 POINT (-100 59)
#> 2  POINT (-90 75)
#> 3 POINT (-120 62)

sf::st_join(my_x, my_y)
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -129 ymin: 56 xmax: -80 ymax: 69
#> Geodetic CRS:  WGS 84
#>          geometry
#> 1 POINT (-129 56)
#> 2  POINT (-80 69)
#> 3 POINT (-100 60)

# Doesn't work. it just does all combinations...
sf::st_union(my_x, my_y)
#> Simple feature collection with 9 features and 0 fields
#> Geometry type: MULTIPOINT
#> Dimension:     XY
#> Bounding box:  xmin: -129 ymin: 56 xmax: -80 ymax: 75
#> Geodetic CRS:  WGS 84
#>                         geometry
#> 1 MULTIPOINT ((-100 59), (-12...
#> 2 MULTIPOINT ((-80 69), (-100...
#> 3 MULTIPOINT ((-100 60), (-10...
#> 4 MULTIPOINT ((-90 75), (-129...
#> 5 MULTIPOINT ((-80 69), (-90 ...
#> 6 MULTIPOINT ((-90 75), (-100...
#> 7 MULTIPOINT ((-120 62), (-12...
#> 8 MULTIPOINT ((-80 69), (-120...
#> 9 MULTIPOINT ((-100 60), (-12...


# I would like a linestring sf object of length 3

Created on 2023-11-09 with reprex v2.0.2

Maybe

st_union(my_x, my_y, by_feature= TRUE)

could work?

Similar to sf::st_distance()

Thanks

[edzer]

I would suggest to do this with by_feature = TRUE; right now by_feature is ignored if y is not missing.

[olivroy]

Sounds like a good idea. Thanks for the clarification. I updated my feature request accordingly to reflect that.

[edzer]

The pairwise geometries you're looking for you can now get with st_union(st_geometry(my_x), st_geometry(my_y)), it's right now pretty inconsistent, at first glance, how st_union() works for sf and for sfc objects.

[olivroy]

It works as you suggest, but I would have never thought to do it this way.

I think this could be added to the st_union docs?

Edit: The st_geometry option has the disadvantage of dropping all non-geometry columns from my_x and my_y. It is pretty easy to recuperate in another operation, but nonetheless important to point out!

For ref

new_line_shp <- st_union(st_geometry(my_x), st_geometry(my_y)) |> 
  sf::st_cast("LINESTRING") |> 
  sf::st_as_sf() |> 
  dplyr::rename(geometry = x)


# to get the geometry column at the end.
bind_cols(
  my_x |> sf::st_drop_geometry(),
  new_line_shp) |> 
  sf::st_as_sf()