如何在R的ggplot2/tmap/mapsf中绘制站点空间对应的数据椭圆

思路：将监测站坐标转为sf空间对象，为每个站点计算污染物数据的椭圆几何（利用ellipse包生成椭圆坐标），最后在地图底图上叠加站点点和椭圆多边形。

步骤与代码 #

1. 加载依赖包

library(ggplot2)
library(sf)
library(dplyr)
library(ellipse)

2. 模拟数据（替换为你的真实数据）

# 监测站基础信息：ID、坐标
stations <- tibble(
  station_id = 1:10,
  lon = runif(10, 116, 117),  # 模拟经度
  lat = runif(10, 39, 40)     # 模拟纬度
)

# 污染物样本数据：每个站点对应多组PM2.5、NO2测量值
pollutants <- stations %>%
  slice(rep(1:n(), each = 50)) %>%  # 每个站点50个样本
  mutate(
    pm25 = rnorm(n(), mean = 30 + station_id*2, sd = 5),
    no2 = rnorm(n(), mean = 20 + station_id*1.5, sd = 4)
  )

3. 计算每个站点的椭圆几何并转为sf对象

# 分组计算椭圆坐标
station_ellipses <- pollutants %>%
  group_by(station_id) %>%
  summarise(
    # 计算污染物的协方差矩阵
    cov_mat = list(cbind(pm25, no2) %>% cov()),
    # 计算污染物均值（椭圆中心）
    mean_pm25 = mean(pm25),
    mean_no2 = mean(no2),
    # 生成椭圆坐标点（置信水平95%）
    ellipse_coords = list(ellipse(cov_mat[[1]], centre = c(mean_pm25, mean_no2), level = 0.95))
  ) %>%
  left_join(stations, by = "station_id") %>%
  # 将椭圆坐标转为sf多边形，同时关联站点地理坐标
  rowwise() %>%
  mutate(
    # 注意：缩放比例需根据实际污染物数值和地图范围调整
    ellipse_geo = list(
      st_polygon(list(
        cbind(
          lon + (ellipse_coords[,1] - mean_pm25)/20,  # 适配经度维度
          lat + (ellipse_coords[,2] - mean_no2)/15    # 适配纬度维度
        )
      ))
    )
  ) %>%
  st_as_sf(sf_column_name = "ellipse_geo")

# 监测站转为sf点对象
stations_sf <- stations %>% st_as_sf(coords = c("lon", "lat"), crs = 4326)

4. 绘制叠加图

ggplot() +
  # 添加地图底图（可替换为你的研究区边界）
  geom_sf(data = st_as_sf(maps::map("state", plot = FALSE, fill = TRUE)), fill = "#f0f0f0", color = "gray") +
  # 叠加椭圆多边形
  geom_sf(data = station_ellipses, fill = "red", alpha = 0.3, color = "darkred") +
  # 叠加监测站点
  geom_sf(data = stations_sf, color = "blue", size = 2) +
  # 自定义样式
  theme_minimal() +
  labs(title = "监测站污染物数据椭圆空间叠加图", x = "经度", y = "纬度")

思路：tmap支持分层绘制空间要素，先加载底图，再依次叠加站点点和预先生成的椭圆多边形，适合快速生成交互式或静态地图。

步骤与代码 #

1. 加载依赖包

library(tmap)
library(sf)
library(dplyr)
library(ellipse)

2. 复用之前的模拟数据和椭圆计算代码（同ggplot2部分的stations、pollutants、station_ellipses、stations_sf）

3. 绘制叠加图

# 设置tmap模式（静态"plot"或交互式"view"）
tmap_mode("plot")

tm_shape(st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))) +
  tm_fill("#f0f0f0") +
  tm_borders("gray") +
# 叠加椭圆
tm_shape(station_ellipses) +
  tm_fill("red", alpha = 0.3) +
  tm_borders("darkred") +
# 叠加站点
tm_shape(stations_sf) +
  tm_dots(col = "blue", size = 0.5) +
# 添加标题和样式
tm_layout(title = "监测站污染物椭圆空间叠加", legend.show = FALSE)

思路：mapsf专注于空间数据可视化，通过mf_map分层绘制底图、椭圆和站点，语法简洁，适合专题地图制作。

步骤与代码 #

1. 加载依赖包

library(mapsf)
library(sf)
library(dplyr)
library(ellipse)

2. 复用模拟数据和椭圆计算代码（同前）

3. 绘制叠加图

# 初始化地图
mf_init(st_as_sf(maps::map("state", plot = FALSE, fill = TRUE)))

# 绘制底图
mf_map(st_as_sf(maps::map("state", plot = FALSE, fill = TRUE)), fill = "#f0f0f0", border = "gray")

# 绘制椭圆
mf_map(station_ellipses, fill = "red", alpha = 0.3, border = "darkred")

# 绘制监测站点
mf_map(stations_sf, pch = 20, col = "blue", cex = 1.5)

# 添加标题
mf_title("监测站污染物数据椭圆空间叠加")

关键注意事项 #

• 椭圆缩放：代码中对污染物数据的缩放比例（如/20、/15）需根据实际污染物数值范围和地图坐标范围调整，确保椭圆大小显示合理。
• 坐标系统：确保所有空间对象使用统一的CRS（示例用WGS84，即EPSG:4326），若使用投影坐标（如UTM），缩放逻辑需对应调整。
• 椭圆参数：可通过ellipse()函数的level参数调整置信水平（默认0.95），或自定义协方差矩阵生成不同类型的椭圆。

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