R sf package: Error in subsetting geometries - "invalid subscript type 'list'"

Question

I'm working with the sf package in R to process a set of point geometries. My objective is to find neighbors for each point within a certain distance. However, I'm encountering an error when attempting to subset geometries.

Here's a simplified version of my code:

library(sf)

# Load the shapefile
points_sf <- st_read("path_to_my_shapefile.shp")

# Define a function to get points within a buffered geometry
get_points_within <- function(geometry, dataset, buffer_distance) {
  buffered_geometry <- st_buffer(geometry, dist = buffer_distance)
  return(dataset[st_within(st_sfc(geometry, crs = st_crs(dataset)), buffered_geometry), ])
}

# Extract geometries from the sf object
geometries <- st_geometry(points_sf)

# Test the function with an individual geometry
single_geometry <- geometries[1]
neighbors <- get_points_within(single_geometry, points_sf, 6)

When running the above code, I get the following error:

Error in `[.default`(geom, i) : invalid subscript type 'list'

From the traceback, it seems that the issue is related to the st_within() function call within the get_points_within function.

Can someone help me understand why this error is occurring and how to resolve it?

Answer 1

st_within() and I think all topological relation functions return a sparse matrix (in a list: see https://r.geocompx.org/spatial-operations.html#topological-relations).

You can have multiple points within your buffer so it is the correct default behavior. In your code you are trying to subset a list (only one dimension) while using two (a[one,two])

As an example:

# from https://r.geocompx.org/spatial-operations.html#topological-relations
polygon_matrix = cbind(
  x = c(0, 0, 1, 1,   0),
  y = c(0, 1, 1, 0.5, 0)
)
polygon_sfc = st_sfc(st_polygon(list(polygon_matrix)))

# create points
point_df = data.frame(
  x = c(0.2, 0.7, 0.4),
  y = c(0.1, 0.2, 0.8)
)
point_sf = st_as_sf(point_df, coords = c("x", "y"))

# I reversed polygon_sfc and point_sf to have a quick working example

st_intersects(polygon_sfc, point_sf)
#Sparse geometry binary predicate
#list of length 1, where the
#predicate was `intersects'
# 1: 1, 3

bob <- st_intersects(polygon_sfc, point_sf)
bob[1,]
# Error in bob[1, ] : incorrect number of dimensions
# the trouble you are runing into

If you just need to know that at least one point is in your buffer you can just add sparse = FALSE

st_intersects(polygon_sfc, point_sf, sparse = FALSE)
#   [,1]  [,2] [,3]
#[1,] TRUE FALSE TRUE
point_df[st_intersects(polygon_sfc, point_sf, sparse = FALSE) , ]