st_intersection to spatialpolygon dataframe not working

Question

I am trying to get the intersection of two shapefiles (census tracts that fall within the boundaries of certain metropolitan areas). I am able to successfully get the intersecting features, however when I try to convert the output of sf_intersect to a SpatialPolygonsDataframe I get the error:

"Error in as_Spatial(from) : conversion from feature type sfc_GEOMETRY to sp is not supported"

This is my code:

library(sf)
library(dplyr)
library(tigris)
library(sp)

#download shapefiles corresponding to metro areas 
metro_shapefiles<-core_based_statistical_areas(cb = FALSE, year = 2016)
#convert to sf and filter
metro_shapefiles<-st_as_sf(metro_shapefiles)%>%filter(GEOID==31080 )
#Data for California
census_tracts_california<-tracts(state="CA",year=2016)
census_tracts_california<-st_as_sf(census_tracts_california)

#INTERSECT AND CONVERT BACK TO SP
census_tracts_intersected1<-st_intersection(census_tracts_california,
                                            metro_shapefiles)

#back to spatial
census_tracts_intersected1<-as(census_tracts_intersected1,"Spatial")

Answer 1

The error message is telling you you can't convert an sfc_GEOMETRY to a Spatial object. There is no sp equivalent object.

In your intersection result you have a mixture of geometries (and hence, you're returned an sfc_GEOMETRY as your 'geometry'). You can see all the geometries here:

types <- vapply(sf::st_geometry(census_tracts_intersected1), function(x) {
    class(x)[2]
}, "")

unique(types)
# [1] "POLYGON"         "MULTILINESTRING" "MULTIPOLYGON"

If you want, you can extract each type of geometry, and convert those to SP individually:

lines <- census_tracts_intersected1[ grepl("*LINE", types), ]
polys <- census_tracts_intersected1[ grepl("*POLYGON", types), ]

spLines <- as(lines, "Spatial")
spPolys <- as(polys, "Spatial")

---

Additional Information

I mentioned in the comments you could use st_join. However, this may not give you the result you want. Within sf library there are the geometric binary predicates, such as ?st_intersects, and geometric operations such as ?st_intersection

The predicates return a sparse (default) or dense matrix telling you with which geometry of y each geometry of x intersects. If you use this within st_join, it will return the (original) geometries that intersect, rather than the sparse matrix.

Whereas the operations (e.g. st_intersection) will compute the intersection, and return new geometries.

Example use

The predicates (st_intersects) can be used inside st_join, and they will return the original geometries which 'intersect'

sf_join <- sf::st_join(census_tracts_california, metro_shapefiles, join = st_intersects)

In this case this gives a single type of object

types <- vapply(sf::st_geometry(sf_join), function(x) {
  class(x)[2]
}, "")

unique(types)
# [1] "MULTIPOLYGON"

## so you can convert to a Spatial object
spPoly <- as(sf_join, "Spatial")

But you need to decide if the result of st_intersect is what you're after, or whether you need the new geometries given by st_intersection.

Further reading

• information on each join is on the sf blog.

• spatial predicates and examples of what the different operations do are on wikipedia (with some good illustrations)

---

Credit to user @lbussett for their description on the difference between st_intersect and st_intersection

Answer 2

Conversion to a Spatial object can't handle mixed lines and polygons. After the intersection you can extract just the polygons (and discard any lines) using:
st_collection_extract("POLYGON")

Your example doesn't seem to fail anymore, so I've created a new example that intersects two polygons, with a shared side. This results in an intersection output of a polygon and a line. On the second attempt I've piped the intersection through the st_collection_extract() function prior to successful conversion to a Spatial object.

library(sf)
library(dplyr)
library(sp)

#Create some boxes
BoxA <- st_polygon(list(cbind(c(0,0,2,2,0),c(0,2,2,0,0))))
BoxB <- st_polygon(list(cbind(c(1,1,3,3,1),c(1,3,3,1,1))))
BoxC <- st_polygon(list(cbind(c(2,2,4,4,2),c(0,2,2,0,0))))

#Create a funny shaped union to help demonstrate the intersection issue
BoxAB <- st_union(BoxA,BoxB)

plot(BoxAB)
plot(BoxC,add=TRUE,border="blue")

Example polygons to intersect

#Intersect of BoxAB with BoxC results in a line and a polygon
BoxIntersects<-st_intersection(BoxAB,BoxC)
plot(BoxIntersects)

Intersection made up of a polygon and a line

#back to spatial fails
SpatialVersionOfIntersects<-as(BoxIntersects,"Spatial")

Error in .as_Spatial(from, cast, IDs) : conversion from feature type sfc_GEOMETRY to sp is not supported

#Intersect again, but this time extract only the polygons
BoxIntersects<-st_intersection(BoxAB,BoxC) %>% st_collection_extract("POLYGON")

#back to spatial suceeds!
SpatialVersionOfIntersects<-as(BoxIntersects,"Spatial")