Why does this R code that maps points into polygons cause a segmentation fault?

Question

As in this question/answer, I have a data frame of points (which also includes some necessary, non-coordinate data columns) and a shapefile. In this case the shapes are US states. The last line of the following code causes a segmentation fault:

library(dplyr)
library(ggplot2)
library(magrittr)
library(sf)
library(USAboundaries)

# Load state shapefile from 1990 Census
map <- USAboundaries::us_states(map_date = "1991-01-01", resolution = "high", 
                                states = c("CA", "OR", "WA", "NV")) %>%
  select(state_name, geometry)

map %>%
  ggplot() +
  geom_sf(aes(fill = state_name))


# Firm/business identifiers, coordinates, and book values
firms <- data.frame(
  "firm_id" = c("A", "B", "C"),
  "lon" = c(41.245696, 38.573965, 44.235577),
  "lat" = c(-118.484406, -121.264677, -115.873234),
  "book_value" = c(15174, 482020, 238592)
)

pnts_sf <- do.call("st_sfc",c(lapply(1:nrow(firms), 
                                     function(i) {st_point(as.numeric(firms[i, c("lon", "lat")]))}), 
                              list("crs" = 4326)))

pnts_trans <- st_transform(pnts_sf, 2163)
map_trans <- st_transform(map, 2163)

# intersect and extract state name
firms$state <- apply(st_intersects(map_trans, pnts_trans, sparse = FALSE), 2, 
                     function(col) { 
                       map_trans[which(col), ]$state_name
                     })

The traceback is straightforward:

Traceback:
 1: CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern,     prepared)
 2: st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared)
 3: st_intersects(map_trans, pnts_trans, sparse = FALSE)
 4: apply(st_intersects(map_trans, pnts_trans, sparse = FALSE), 2,     function(col) {        map_trans[which(col), ]$state_name    })
An irrecoverable exception occurred. R is aborting now ...
Segmentation fault

I'm using R 3.5.2 and sf 0.7-3, if it helps. I'm assuming this is a bug in my code and not the libraries, but can someone replicate this? Is there something wrong with my code?

Answer 1

Running this I also got the same error. Usually seg faults are bugs in underlying C code from packages. Valid R code shouldn't cause a seg fault, even if it is wrong.

The st_intersects was at fault here so I looked at the objects being passed to it:

> pnts_trans
Geometry set for 3 features  (with 3 geometries empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    2163
proj4string:    +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
POINT EMPTY
POINT EMPTY
POINT EMPTY

Empty points? That's odd. Before the transformation pnts_sf was valid points. Hmmm.... Wait a minute...

"lon" = c(41.245696, 38.573965, 44.235577),
"lat" = c(-118.484406, -121.264677, -115.873234),

You can't have a latitude less than -90 (that's the south pole). Looks like you've mixed latitude and longitude (not the first time anyone has done that...).

With those latitudes st_transform silently converts the points to an empty geometry, and then st_intersect with empty point geomtries causes the seg fault. So its a combination of your mistake plus a bug in sf. You should file a bug report for sf when using st_intersect with EMPTY point geometries.