sp × tidyverse test: rbind

This is a set of tests related to a Github issue for the sf package.

library(tidyverse)
library(sf)

Make test objects

Make a polygon, a multipolygon, and two dataframes that can be combined with the geometries to create simple feature objects.

# Polygon geometry: pol
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
# 
# Multi-polygon geometry: mpol
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5)))
# 
# sf object w/ polygon: sf_pol
df1 <- 
        data.frame(VAR1 = 'foo',VAR2 = 'bar') %>% 
        mutate(geom = st_sfc(pol))
sf_pol <- st_as_sf(df1) 
# 
# sf object w/ multi-polygon: sf_mpol
df2 <- 
        data.frame(VAR1 = 'foo',VAR2 = 'bar') %>% 
        mutate(geom = st_sfc(mpol))
sf_mpol <- st_as_sf(df2)

Test 1: rbind two sf objects with different geometry classes

# Check the sf object metadata: sf_pol
sf_pol

Simple feature collection with 1 feature and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 4
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR2                           geom
1  foo  bar POLYGON((0 0, 1 0, 3 2, 2 4...

# Check the sf object metadata: sf_mpol
sf_mpol

Simple feature collection with 1 feature and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 4
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR2                           geom
1  foo  bar MULTIPOLYGON(((0 0, 1 0, 3 ...

# Attempt to combine the object with `rbind()`
sf_both <- rbind(sf_pol,sf_mpol)
sf_both

Simple feature collection with 2 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 4
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR2                           geom
1  foo  bar POLYGON((0 0, 1 0, 3 2, 2 4...
2  foo  bar MULTIPOLYGON(((0 0, 1 0, 3 ...

Looking at the metadata above, I’m surprised to see geometry type: POLYGON. Based on the sf package vignette, I would have guessed that the geometry of this new simple feature would be of the GEOMETRY class (but perhaps that was wishful thinking).

Let’s see what happens when we plot this object:

try(plot(sf_both))

Error in eval(substitute(expr), envir, enclos) : 
  not compatible with requested type

It seems that plot() cannot handle this mixed-type sfc. It would be ideal if rbind() could detect the presence of different types of sfg’s and coerce them into a GEOMETRY sfc.

Otherwise, users need to find work-arounds. The following is one way to manually build the GEOMETRY sfc by looping over the two sfg’s:

sfc_new <- st_sfc()
        for(i in seq_along(st_geometry(sf_pol))){
                 
                 geom <- st_geometry(sf_pol)[[i]]
                 sfc_new[[i]] <- geom
        }
        for(i in seq_along(st_geometry(sf_mpol))){
                 
                 geom <- st_geometry(sf_mpol)[[i]]
                 sfc_new[[i + nrow(sf_pol)]] <- geom
                 
        }
sf_both_new <- 
        rbind(sf_pol,
              sf_mpol) %>% 
        mutate(geom = sfc_new)
sf_both_new

Simple feature collection with 2 features and 2 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR2                           geom
1  foo  bar POLYGON((0 0, 1 0, 3 2, 2 4...
2  foo  bar MULTIPOLYGON(((0 0, 1 0, 3 ...

try(plot(sf_both_new))

So that work-around returns the desired results, but it wouldn’t scale very well to accomodate the merger of many sf objects and it certainly breaks the continuity and readability of a %>% workflow (see princples 2 and 4 of the Tidy Tools Manifesto).

Test 2: bind sf objects with non-matching columns and different geometry types

# sf object w/ polygon: sf_pol2
df3 <- data.frame(VAR1 = 'foo',VAR2 = 'bar') 
df3$geom <- st_sfc(pol)
sf_pol2 <- st_as_sf(df1) 
# 
# sf object w/ multi-polygon: sf_mpol2 (new column: VAR3)
df4 <- data.frame(VAR1 = 'foo',VAR3 = 'bar') 
df4$geom <- st_sfc(mpol)
sf_mpol2 <- st_as_sf(df4)
# 
# Try rbind()
sf_pol2;sf_mpol2

Simple feature collection with 1 feature and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 4
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR2                           geom
1  foo  bar POLYGON((0 0, 1 0, 3 2, 2 4...
Simple feature collection with 1 feature and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 4
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR3                           geom
1  foo  bar MULTIPOLYGON(((0 0, 1 0, 3 ...

try(rbind(sf_pol2,sf_mpol2))

Error in match.names(clabs, names(xi)) : 
  names do not match previous names

# 
# Try dplyr::bind_rows()
# NOTE: this throws an error, even with the help of try() 
# 
# try(dplyr::bind_rows(sf_pol2,sf_mpol2))
# 
# Work-around
df_both_new2 <- dplyr::bind_rows(sf_pol2[,1:2],
                                 sf_mpol2[,1:2])
df_both_new2$geom <- sfc_new
sf_both_new2 <- st_sf(df_both_new2)
sf_both_new2; try(plot(sf_both_new2))

Simple feature collection with 2 features and 3 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    NA
proj4string:    NA
  VAR1 VAR2 VAR3                           geom
1  foo  bar <NA> POLYGON((0 0, 1 0, 3 2, 2 4...
2  foo <NA>  bar MULTIPOLYGON(((0 0, 1 0, 3 ...

dplyr::bind_rows() fails because it doesn’t recognize the sf class. However, stripping the sf objects of their sfc’s, passing them to bind_rows(), and then adding the manually-built GEOMETRY sfc returns the desired result.

Thoughts

Perhaps these tests make the case for the following enhancements to sf:

1. an sf-compatible version of dplyr::bind_rows() (provided it’s even possible - see Hadley’s comment)
2. an sf function for combining sfc’s with different geometry types into a sfc with GEOMETRY type sfg’s