如何裁剪关联空间数据？解决sf与data.frame兼容报错问题

问题背景 #

现有两个数据集：

• df1：包含纬度（Latitude）和经度（Longitude）的data.frame，数据结构如下：

df1 = structure(list(Latitude = c(44.11, 45.78, 44.49, 46.87, 44.81, 
44.11, 44.45, 45.78, 44.79, 45.22, 47.3, 44.24, 44.33, 44.96, 
44.98, 44.73, 44.36, 45.14, 45.98, 46.59, 46.53, 44.48, 45.09, 
44.31, 45.79, 44.74, 45.02, 44.71, 44.55, 46.91, 44.58, 44.1, 
45.39, 44.84, 45.08, 46.16, 44.49, 44.31, 44.37, 44.36, 44.38, 
46.88, 46.21, 44.31, 44.7, 44.68), Longitude = c(-70.15, -87.08, 
-73.12, -68.03, -68.8, -70.15, -71.17, -87.08, -85.64, -67.26, 
-68.57, -73.49, -71.75, -72.16, -74.82, -73.44, -74.11, -87.6, 
-86.21, -87.4, -84.33, -83.33, -83.41, -85.39, -84.71, -84.69, 
-84.66, -85.64, -87.92, -67.99, -70.54, -70.17, -68.52, -68.72, 
-69.86, -67.82, -73.2, -69.68, -69.75, -69.74, -69.79, -68.02, 
-67.79, -69.68, -67.57, -67.64)), row.names = c(NA, -46L), class = "data.frame")

• usa：sf格式的美国州界多边形数据，生成代码如下：

library(sf)
library(rnaturalearth)

world <- ne_countries(scale = "medium", returnclass = "sf")
usa = filter(world,admin =="United States of America")
usa <- st_as_sf(maps::map("state", fill=TRUE, plot =FALSE))

尝试用SpatialPointsDataFrame筛选境内点时出现如下报错：

Error in h(simpleError(msg, call)) :
error in evaluating the argument 'y' in selecting a method for function 'over': internal problem in as(): “sf” is(object, "data.frame") is TRUE, but the metadata asserts that the 'is' relation is FALSE"

解决方案 #

报错根源是混用了sp包（SpatialPointsDataFrame）和sf包的对象，导致类型冲突。全程使用sf包的工具即可避免此类问题，步骤如下：

1. 修正美国边界数据的坐标系 #

首先确保usa的坐标系与df1的经纬度坐标系（WGS84，EPSG:4326）一致：

library(sf)
library(dplyr)
library(maps)

# 重新生成并规范美国州界的sf对象
usa <- st_as_sf(maps::map("state", fill = TRUE, plot = FALSE)) %>%
  st_set_crs(4326)  # 设置为WGS84坐标系

2. 将df1转换为sf格式的点数据 #

把普通的data.frame转为sf点对象，指定经纬度列和坐标系：

df1_sf <- st_as_sf(df1, coords = c("Longitude", "Latitude"), crs = 4326)

3. 筛选境内点 #

有两种高效方法：

方法一：用st_within判断归属

生成逻辑向量标记每个点是否在美国境内，再筛选：

# 判断每个点是否位于美国多边形内，sparse=FALSE返回矩阵，取第一列得到逻辑向量
inside_usa <- st_within(df1_sf, usa, sparse = FALSE)[, 1]

# 筛选出美国境内的点
df1_usa <- df1_sf[inside_usa, ]

方法二：用st_intersection直接裁剪

直接获取点与美国边界的交集，结果即为境内点：

df1_usa <- st_intersection(df1_sf, usa)

处理完成后，df1_usa就是位于美国境内的点数据，保留原df1的所有属性，同时带有sf空间属性。

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