如何基于经纬度空间数据绘制Jaccard Index与空间距离对比图

1. 环境准备与数据构造 #

首先安装加载所需依赖包：

# 首次使用先安装
install.packages(c("geosphere", "ggplot2", "reshape2"))
# 加载包
library(geosphere)
library(ggplot2)
library(reshape2)

先准备好符合要求的输入数据：

• 经纬度数据：两列数据框/矩阵，第一列为经度，第二列为纬度，每行对应一个点位
• Jaccard Index矩阵：和点位数量匹配的n*n对称矩阵，对角线为1，值为对应两个点位的Jaccard指数

示例构造测试数据：

set.seed(123)
# 构造10个随机点位的经纬度，范围是东经100-120°，北纬30-40°
lon_lat <- data.frame(
  lon = runif(10, 100, 120),
  lat = runif(10, 30, 40)
)
# 构造对称的10*10 Jaccard矩阵
jaccard_mat <- matrix(runif(100, 0, 1), nrow = 10)
jaccard_mat[lower.tri(jaccard_mat)] <- t(jaccard_mat)[lower.tri(jaccard_mat)]
diag(jaccard_mat) <- 1

2. 用geosphere计算空间距离矩阵 #

geosphere提供多种距离计算函数，常用的两个：

• distHaversine：半正矢公式，假设地球为正球体，计算速度快，适合中小范围点位距离计算
• distVincentyEllipsoid：文森特公式，基于椭球模型，精度更高，适合大范围高精度需求

默认返回结果单位为米，可按需转换为千米：

# 计算成对空间距离矩阵
dist_mat <- distHaversine(lon_lat)
# 单位转换为千米
dist_mat <- dist_mat / 1000

注意：所有geosphere的距离函数都要求输入经纬度顺序为经度在前、纬度在后，顺序颠倒会导致结果完全错误

3. 矩阵匹配与对比绘图 #

两个矩阵都是对称矩阵，我们仅提取下三角非对角线的非重复值配对，避免数据冗余：

# 提取下三角非对角线值
dist_vals <- dist_mat[lower.tri(dist_mat, diag = FALSE)]
jaccard_vals <- jaccard_mat[lower.tri(jaccard_mat, diag = FALSE)]
# 构造绘图数据集
plot_df <- data.frame(
  spatial_distance_km = dist_vals,
  jaccard_index = jaccard_vals
)

绘制散点对比图，可额外添加拟合线查看趋势：

ggplot(plot_df, aes(x = spatial_distance_km, y = jaccard_index)) +
  geom_point(alpha = 0.6, size = 2) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(x = "空间距离(km)", y = "Jaccard指数") +
  theme_bw()

注意事项 #

• 经纬度如果是度分秒格式，需要先转换为十进制度再传入计算函数
• 距离矩阵和Jaccard矩阵的点位顺序必须完全对应，建议计算前核对行名/列名是否一致
• 点位数量超过1000时可抽样部分值绘图，避免散点过度重叠影响可读性

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