使用sf包st_intersects测试点与几何相交遇错求助

问题背景 #

我正尝试使用st_intersects()函数测试自有点数据是否落在地图数据的几何图形范围内。

使用的数据：

• 第一组：加州消防局野火区域shapefile
• 第二组：另一组加州消防局野火区域shapefile

当前先用以下测试点数据验证是否落在上述shapefile的多边形中：

la_test_points <- data.frame(y = runif(1000, 33.6, 34.8), x = runif(1000, -119, -117.6))

将地图数据与点数据叠加后显示正常，但执行相交测试时出现问题。

出错代码及错误信息 #

使用st_intersects的错误 #

# 转换shapefile地图数据的坐标系
la_fire_sra <- st_transform(st_as_sf(la_fire_sra), crs = 3857)
# 将测试点转为sf对象并匹配坐标系
la_test_points_merged <- st_as_sf(la_test_points, coords = c('y', 'x'), crs = st_crs(la_fire_sra))
# 测试点是否落在shapefile的任一几何图形内
la_test_points_merged <- la_test_points_merged %>%
  mutate(intersection = st_intersects(geometry, la_fire_sra))

执行后打印la_test_points_merged时出现错误：

> la_test_points_merged
Simple feature collection with 1000 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 33.60155 ymin: -118.9959 xmax: 34.79907 ymax: -117.6015
Projected CRS: WGS 84 / Pseudo-Mercator
First 10 features:
Error in xj[i, , drop = FALSE] : incorrect number of dimensions

使用st_intersection的错误 #

改用st_intersection()函数时出现另一个错误：

> la_test_points_merged <- la_test_points_merged %>%
+   mutate(intersection = st_intersection(geometry, la_fire_sra))
Error in `stopifnot()`:
! Problem while computing `intersection = st_intersection(geometry, la_fire_sra)`.
x `intersection` must be size 1000 or 1, not 0.
Run `rlang::last_error()` to see where the error occurred.

我期望得到的结果是判断每个测试点是否被la_fire_sra的任一几何图形包含，即常规的点-in-polygon判断结果。请问如何修改代码解决该问题？

解决方案 #

错误原因分析 #

1. 坐标顺序错误：创建点sf对象时，coords = c('y', 'x')搞反了地理坐标的顺序（经度对应x，纬度对应y），导致点的位置完全偏离多边形区域，最终无相交结果。
2. st_intersects使用方式不当：默认返回的是列表格式的匹配索引，直接放入mutate会导致数据结构不兼容；而st_intersection是返回相交的几何对象，并非布尔判断，无相交时会返回空值触发报错。

修改后的代码 #

# 1. 读取并转换多边形数据坐标系
la_fire_sra <- st_as_sf(la_fire_sra) %>% st_transform(crs = 3857)

# 2. 创建测试点sf对象：修正坐标顺序，先x（经度）后y（纬度）
la_test_points <- data.frame(
  x = runif(1000, -119, -117.6),  # 经度
  y = runif(1000, 33.6, 34.8)     # 纬度
)
# 先指定原始WGS84坐标系，再转换到与多边形一致的3857
la_test_points_merged <- st_as_sf(la_test_points, coords = c('x', 'y'), crs = 4326) %>% 
  st_transform(crs = st_crs(la_fire_sra))

# 3. 方式一：用lengths判断是否有匹配的多边形，得到布尔值
la_test_points_merged <- la_test_points_merged %>%
  mutate(is_inside = lengths(st_intersects(geometry, la_fire_sra)) > 0)

# 方式二：用sparse=FALSE返回逻辑矩阵，取行逻辑或得到布尔值
la_test_points_merged <- la_test_points_merged %>%
  mutate(is_inside = rowSums(st_intersects(geometry, la_fire_sra, sparse = FALSE)) > 0)

# 查看结果
head(la_test_points_merged)

关键说明 #

• 坐标顺序修正：地理坐标必须遵循「经度（x）在前，纬度（y）在后」的顺序，否则点的空间位置会完全错误，导致与多边形无交集。
• st_intersects正确用法：
  • 默认sparse=TRUE返回列表，每个元素是当前点匹配的多边形索引，通过lengths()判断是否有匹配即可得到「是否在多边形内」的布尔值。
  • sparse=FALSE返回逻辑矩阵，每行对应一个点，每列对应一个多边形，rowSums(...)>0表示该点至少落在一个多边形内。
• st_intersection的适用场景：该函数用于提取相交的几何对象，而非判断是否相交，不适合用来生成布尔型的点-in-polygon结果。

内容的提问来源于stack exchange，提问作者qwerty