如何判断R中SpatialPointsDataFrame的点是否落在指定50km缓冲区内

你现有代码的核心问题是WGS84地理坐标系的单位是十进制度，直接给gBuffer传入米为单位的宽度会出现单位不匹配的问题，才会生成异常大的缓冲区。另外不需要额外生成缓冲区文件，直接计算两组坐标的距离就能实现你的需求，操作更简单效率更高。

优化实现方案（基于sf包，sp/rgdal/rgeos已停止维护，sf为目前R空间生态主流工具） #

# 加载所需包
library(dplyr)
library(sf)

# 1. 读取并清理空值数据
FDI <- read.csv("FDI_harmonized.csv") %>%
  filter(!is.na(oc_lat), !is.na(oc_lng), !is.na(oc_lat_st), !is.na(oc_lng_st))

# 2. 分别构造两组坐标的空间点对象，指定WGS84坐标系（EPSG:4326和你原有proj参数完全一致）
# 第一组：oc_lng/oc_lat对应点
pts_oc <- st_as_sf(FDI, coords = c("oc_lng", "oc_lat"), crs = 4326)
# 第二组：oc_lng_st/oc_lat_st对应点
pts_oc_st <- st_as_sf(FDI, coords = c("oc_lng_st", "oc_lat_st"), crs = 4326)

# 3. 逐行计算两组点的大圆距离（默认单位为米），判断是否在50km范围内
FDI$within_50km <- as.numeric(st_distance(pts_oc, pts_oc_st, by_element = TRUE)) <= 50 * 1000

说明 #

• by_element = TRUE参数会按行逐对计算两个点集的距离，不会生成全量距离矩阵，计算效率更高
• 最终生成的within_50km列就是你需要的布尔标记：TRUE代表第二组坐标在第一组坐标的50km范围内，FALSE反之，如需转成1/0标记，直接调用as.integer(FDI$within_50km)即可

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