如何在R中通过坐标范围筛选几何列获取特定区域河流数据

问题重现 #

我正在制作某区域的河流地图，从HydroRivers获取了河流数据，想通过指定坐标范围子集化数据，但以下代码运行报错：

library(sf)
library(tidyverse)

river_data <- read_sf("~/HydroRIVERS_v10_eu.shp")
x_range <- c(25.7, 45.3)
y_range <- c(35.7, 42.3)

my_selection <- river_data$geometry[st_intersects(river_data$geometry, st_bbox(c(x_range, y_range))),]

数据列信息：

[1] "HYRIV_ID"   "NEXT_DOWN"  "MAIN_RIV"   "LENGTH_KM"  "DIST_DN_KM" "DIST_UP_KM"
[7] "CATCH_SKM"  "UPLAND_SKM" "ENDORHEIC"  "DIS_AV_CMS" "ORD_STRA"   "ORD_CLAS"  
[13] "ORD_FLOW"   "HYBAS_L12"  "geometry"  

Geometry列详情：

Geometry set for 6 features 
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 56.49375 ymin: 81.33125 xmax: 59.30625 ymax: 81.78958
Geodetic CRS:  WGS 84
First 5 geometries:
LINESTRING (59.2625 81.78958, 59.30625 81.78958)
LINESTRING (58.24167 81.77708, 58.13542 81.7770...
LINESTRING (57.50833 81.51042, 57.44375 81.5104...
LINESTRING (57.07917 81.50417, 57.06875 81.5145...
LINESTRING (56.55208 81.33125, 56.54792 81.3354...

错误原因 #

1. st_bbox()生成的是边界框对象，而非sf几何对象，st_intersects()无法直接处理该类型输入
2. 直接索引river_data$geometry会丢失原数据的属性列，不符合保留目标区域完整行数据的需求

正确解决方案 #

方法1：使用st_crop()快速裁剪（推荐） #

st_crop()是sf包专门用于按边界框裁剪空间数据的函数，操作简单且保留所有属性列：

library(sf)
library(tidyverse)

river_data <- read_sf("~/HydroRIVERS_v10_eu.shp")
# 将坐标范围转换为可用于裁剪的sf几何对象
crop_extent <- st_bbox(c(xmin = 25.7, ymin = 35.7, xmax = 45.3, ymax = 42.3)) %>% 
  st_as_sfc()

# 裁剪数据，保留目标区域内及相交的河流要素
cropped_rivers <- st_crop(river_data, crop_extent)

方法2：用st_filter()筛选相交要素 #

适合需要自定义空间筛选逻辑的场景，比如精确控制保留部分相交/完全包含的要素：

# 构建目标区域的多边形边界框
target_bbox <- st_bbox(c(xmin = 25.7, ymin = 35.7, xmax = 45.3, ymax = 42.3)) %>% 
  st_as_sfc()

# 筛选与边界框相交的河流要素
selected_rivers <- river_data %>% 
  st_filter(target_bbox, .predicate = st_intersects)

方法3：直接按坐标范围过滤（快速但有局限） #

仅筛选要素自身边界框完全落在目标范围内的河流，会漏掉跨边界的部分线要素，适合快速初步过滤：

# 提取每个要素的边界框信息
river_bbox_df <- st_bbox(river_data) %>% as.data.frame()

# 按坐标范围筛选
filtered_rivers <- river_data %>% 
  filter(
    xmin >= 25.7,
    xmax <= 45.3,
    ymin >= 35.7,
    ymax <= 42.3
  )

注意事项 #

• 确保目标区域的坐标参考系（CRS）与river_data的CRS一致（示例中为WGS 84），若不一致需先用st_transform()转换CRS
• 方法1和方法2会保留所有与目标区域相交的线要素，更符合地图制作的实际需求

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