如何用terra包将含Quad Key的Data Frame转为栅格（弃用quadkeyr）

我有一个包含Quad Key（四叉键）和活动指数相关值的Data Frame，希望将其转换为栅格。由于quadkeyr包处理加拿大全区域逐小时规模的数据时速度过慢，因此仅打算使用terra包实现转换。以下是数据样例：

structure(list(geography = c("021302200001103223", "021302200001123321", 
"021302200001111311", "021302200001110301", "021302200001100333", 
"021302200001111211"), xlon = c(-112.35512, -112.34962, -112.32491, 
-112.33864, -112.35786, -112.3304), xlat = c(50.72298, 50.70994, 
50.73254, 50.73254, 50.72994, 50.73254), bounds = c("-112.35580,50.72255,-112.35443,50.72342", 
"-112.35031,50.70950,-112.34894,50.71037", "-112.32559,50.73211,-112.32422,50.73298", 
"-112.33932,50.73211,-112.33795,50.73298", "-112.35855,50.72950,-112.35718,50.73037", 
"-112.33109,50.73211,-112.32971,50.73298"), start_date = c("2021-07-01", 
"2021-07-01", "2021-07-01", "2021-07-01", "2021-07-01", "2021-07-01" 
), end_date = c("2021-07-31", "2021-07-31", "2021-07-31", "2021-07-31", 
"2021-07-31", "2021-07-31"), agg_day_period = c(0L, 0L, 1L, 1L, 
1L, 1L), agg_time_period = c(15L, 17L, 10L, 11L, 12L, 15L), activity_index_total = c(0.010026, 
0.034617, 0.02384, 0.051245, 0.023967, 0.027478)), row.names = c(NA, 
6L), class = "data.frame")

我没有原始分辨率的数据，但Data Frame中包含每个要素的bounds信息，请问该如何在不使用专用Quad Key包的前提下完成栅格转换？

解决方案 #

步骤1：解析bounds为空间多边形 #

先把bounds列的字符串格式转换为terra可识别的多边形向量，同时关联数据框中的属性字段：

library(terra)

# 加载数据框（替换为你的实际数据）
df <- structure(...)

# 拆分bounds字符串为数值矩阵，格式为[minlon, minlat, maxlon, maxlat]
bounds_mat <- matrix(as.numeric(unlist(strsplit(df$bounds, ","))), ncol=4, byrow=TRUE)

# 逐个创建多边形并转换为SpatVector
polygon_list <- lapply(1:nrow(bounds_mat), function(i) {
  # 构建多边形闭合坐标：左下→右下→右上→左上→左下
  coords <- matrix(c(
    bounds_mat[i,1], bounds_mat[i,2],
    bounds_mat[i,3], bounds_mat[i,2],
    bounds_mat[i,3], bounds_mat[i,4],
    bounds_mat[i,1], bounds_mat[i,4],
    bounds_mat[i,1], bounds_mat[i,2]
  ), ncol=2, byrow=TRUE)
  polygons(coords)
})
sv <- vect(polygon_list, crs = "EPSG:4326") # 假设坐标系为WGS84，按需调整

# 关联活动指数、时间周期等属性到向量
sv <- cbind(sv, df[, c("activity_index_total", "agg_day_period", "agg_time_period")])

步骤2：创建匹配的栅格模板 #

根据所有多边形的范围和Quad Key单元格大小创建栅格模板：

# 计算单个Quad Key单元格的分辨率（假设所有Quad Key级别一致，取第一个单元格的尺寸）
cell_width <- bounds_mat[1,3] - bounds_mat[1,1]
cell_height <- bounds_mat[1,4] - bounds_mat[1,2]

# 获取所有多边形的整体范围
total_ext <- ext(sv)

# 创建栅格模板
r_template <- rast(total_ext, resolution = c(cell_width, cell_height), crs = crs(sv))

步骤3：多边形属性栅格化 #

使用rasterize函数将活动指数值转换为栅格，支持单波段或多波段（按时间分组）：

# 单波段栅格化：将所有活动指数值转换为一个栅格
r_activity <- rasterize(sv, r_template, field = "activity_index_total", fun = "first")

# 多波段栅格化：按天/小时周期分组生成多波段栅格
time_groups <- unique(df[, c("agg_day_period", "agg_time_period")])
raster_collection <- list()

for(i in 1:nrow(time_groups)) {
  # 筛选对应时间组的多边形
  group_sv <- sv[sv$agg_day_period == time_groups[i,1] & sv$agg_time_period == time_groups[i,2], ]
  # 栅格化当前组
  group_raster <- rasterize(group_sv, r_template, field = "activity_index_total", fun = "first")
  names(group_raster) <- paste0("day_", time_groups[i,1], "_hour_", time_groups[i,2])
  raster_collection[[i]] <- group_raster
}

# 合并为多波段栅格
multi_band_raster <- rast(raster_collection)

关键注意事项 #

• 坐标系匹配：如果你的bounds不是WGS84（EPSG:4326），需要将crs参数改为实际使用的坐标系。
• 分辨率一致性：若数据中存在不同级别的Quad Key，需先按Quad Key长度分组，分别创建对应分辨率的栅格模板再处理。
• 聚合函数：如果存在多边形重叠（同一栅格单元格对应多个Quad Key），可将fun参数改为mean、sum等聚合函数，按需选择。

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