如何在R中计算采样点到值>0.7的人类改造指数栅格单元最短距离？

核心结论 #

完全可以通过转换到UTM坐标系（米单位）来计算距离，最后除以1000即可得到公里单位的结果，UTM坐标系的平面距离计算精度较高，适合这类空间分析需求。

实现步骤与代码 #

1. 加载依赖包

使用sf处理空间矢量数据，raster（或terra）处理栅格数据，dplyr做数据整理：

library(sf)
library(raster)
library(dplyr)

2. 读取并预处理采样点数据

假设你的采样点存储在CSV文件中，包含site_id（站点编号）、lat（纬度）、lon（经度）列：

# 读取采样点并转为空间矢量对象，初始CRS设为WGS84（EPSG:4326）
sample_points <- read.csv("你的采样点文件路径.csv") %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

3. 读取并筛选高人类改造区栅格

读取人类改造指数栅格，提取值>0.7的单元并生成中心点：

# 读取栅格数据（替换为你的栅格文件路径）
hmi_raster <- raster("人类改造指数栅格文件.tif")
# 确保栅格CRS与采样点一致（若未自动识别）
crs(hmi_raster) <- CRS("+init=epsg:4326")

# 提取值>0.7的栅格单元，生成单元中心点的空间对象
high_mod_cells <- rasterToPoints(hmi_raster, fun = function(x) x > 0.7, spatial = TRUE) %>%
  st_as_sf() %>%
  st_centroid()

4. 转换到UTM坐标系

自动匹配采样点所在区域的UTM带（若采样点分散在全球，需分区域处理或改用全球等距坐标系）：

# 获取采样点对应的UTM坐标系
utm_crs <- st_transform(sample_points, crs = "+proj=utm +zone=auto +datum=WGS84 +units=m") %>%
  st_crs()

# 转换采样点和高改造区中心点到UTM
sample_utm <- st_transform(sample_points, crs = utm_crs)
high_mod_utm <- st_transform(high_mod_cells, crs = utm_crs)

5. 计算最短距离并整理结果

计算每个采样点到最近高改造区中心点的距离，转换为公里后合并到原始数据框：

# 计算最短距离（单位：米），转为公里
min_dist_km <- st_nearest_distance(sample_utm, high_mod_utm) %>%
  as.numeric() / 1000

# 合并结果，移除空间属性得到普通数据框
result_df <- sample_points %>%
  st_drop_geometry() %>%
  mutate(到高人类改造区距离_km = min_dist_km)

# 查看结果
head(result_df)

注意事项 #

• 若栅格数据量极大，rasterToPoints可能占用过多内存，推荐使用terra包替代raster，处理效率更高，代码逻辑基本一致。
• 若采样点分布在全球多个UTM带，可按带拆分数据分别计算，再合并结果，避免跨带距离误差。

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