如何在R中生成显示各类别点占比的栅格地图

我们可以用R的terra包（替代旧的raster包，更高效）来实现这个需求，步骤如下：

1. 准备工作 #

首先安装并加载所需包：

# 首次运行安装依赖包
install.packages(c("terra", "dplyr", "viridis"))

# 加载包
library(terra)
library(dplyr)
library(viridis)

2. 数据转换与栅格模板创建 #

将你的点数据转换为矢量对象，并创建匹配范围的栅格模板（可自定义分辨率）：

# 你的示例数据
df <- data.frame(
category = c("A","A","B","A","A","B","A","A","C","B","B","C","D","B","C","B","D","C","B","D"),
X = c(4599984.10552057,4831788.90207249,4840009.6065062,4832038.78315722,4930825.32065239,4416629.65604407,4446435.30204568,4446438.1996941,4446638.22211161,4446635.32438607,4442550.72821205,4442553.62813824,4442750.75801933,4422604.29347903,4422602.42523751,4422805.49222056,4418632.41506483,4418832.47128213,4408657.25277069,4408857.25256069),
Y = c(2881156.36674124,2862062.90168895,3020826.90515788,2864693.62543509,2955978.25574188,3502994.58644907,3502550.9888811,3502750.90387191,3502747.90897273,3502547.99413997,3510607.49600397,3510807.41086243,3510604.49974218,3514904.75866658,3515102.48878262,3515100.46862153,3516963.23902827,3516960.26748617,3519109.88091922,3519310.80154912))

# 转换为SpatVector（terra的矢量格式），这里假设投影是ETRS89 LAEA（欧洲常用投影EPSG:3035），请根据实际数据调整crs参数
points <- vect(df, geom = c("X", "Y"), crs = "EPSG:3035")

# 创建栅格模板：覆盖所有点的范围，设置分辨率（单位与投影一致，这里是米，示例设为10km）
resolution <- 10000  # 可根据需求修改分辨率
r_template <- rast(ext(points), res = resolution, crs = crs(points))

3. 统计栅格单元内的类别占比 #

通过栅格化统计每个单元的总点数和各类别点数，再计算占比：

# 按类别统计每个栅格单元内的点数
category_counts <- rasterize(points, r_template, field = "category", 
                             fun = function(x, ...) table(factor(x, levels = unique(df$category))))

# 将统计结果转为数据框，计算每个类别的占比（处理总点数为0的情况，避免除以0）
ratio_df <- as.data.frame(category_counts, xy = TRUE) %>%
  rename_with(~gsub("V", "", .x), starts_with("V")) %>%
  mutate(total_points = rowSums(across(all_of(unique(df$category))), na.rm = TRUE)) %>%
  mutate(across(all_of(unique(df$category)), ~ifelse(total_points == 0, 0, .x / total_points)))

4. 生成类别占比栅格并可视化 #

将占比数据转回栅格对象，然后分别绘制每个类别的占比地图：

# 生成每个类别的占比栅格
category_ratio_rasts <- lapply(unique(df$category), function(cat) {
  rast(ratio_df, type = "xyz", crs = crs(r_template))[[cat]]
})
names(category_ratio_rasts) <- unique(df$category)

# 绘制所有类别的占比栅格（2x2布局）
par(mfrow = c(2, 2))
for (cat in unique(df$category)) {
  plot(category_ratio_rasts[[cat]], 
       main = paste("类别", cat, "占比"), 
       col = viridis(10), 
       mar = c(2,2,2,2))
  points(points, pch = 16, cex = 0.5, col = "black")  # 叠加原始点作为参考
}
par(mfrow = c(1,1))  # 恢复默认绘图布局

关键说明 #

• 投影匹配：确保点数据和栅格的投影一致，示例用的是欧洲常用的ETRS89 LAEA（EPSG:3035），请根据你的实际数据调整crs参数。
• 分辨率调整：修改resolution参数可以改变栅格单元的大小，单位与投影的坐标单位一致（示例中是米）。
• 空值处理：当栅格单元内没有点时，占比设为0，避免出现NaN。

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