如何用R语言terra包创建不规则栅格？相关数据处理方案咨询

我有一组全球经纬度网格数据集，包含单元格中心的x（经度）、y（纬度）坐标，以及单元格在x、y方向的总长度，采用经纬度CRS。数据示例如下：

head(latlong)
   x      y     x.leng  y.leng   value 
80.221 50.854  72334 16027.8       0.1
80.000 50.214  32132  2404.4       0.2
79.779 50.214  72334 16027.8       0.6
79.664 50.854  77632 23841.8       0.5
79.572 51.548  80556 29873.3       0.3
79.491 52.270  82574 34953.4       0.2

其中x、y为单元格中心坐标，x.leng和y.leng为单元格在对应方向的总长度。由于需要进行基于面积的计算（如提取与多边形重叠的栅格单元格值），无法转换为规则尺寸栅格且不想插值，因此rasterize和rasterXYZ函数不适用。我了解投影后的rast对象可包含不规则尺寸单元格，但不确定如何基于上述数据框创建。是否应考虑创建多要素多边形？还有哪些可行方案？

1. 构建不规则SpatRaster（terra包） #

直接用terra包创建带不规则单元格的栅格对象，完全适配你的面积计算需求：

• 第一步：将米级的单元格长度转换为经纬度度数差（经纬度CRS下距离与纬度相关）：

  library(terra)
  
  # 计算经度方向度数差（1度经度在赤道≈111320米，随纬度升高递减）
  latlong$x_deg <- latlong$x.leng / (111320 * cos(latlong$y * pi/180))
  # 计算纬度方向度数差（1度纬度≈111139米）
  latlong$y_deg <- latlong$y.leng / 111139
  
  # 生成每个单元格的四至边界
  latlong$xmin <- latlong$x - latlong$x_deg/2
  latlong$xmax <- latlong$x + latlong$x_deg/2
  latlong$ymin <- latlong$y - latlong$y_deg/2
  latlong$ymax <- latlong$y + latlong$y_deg/2
  

• 第二步：基于边界创建不规则栅格：

  # 创建SpatRaster对象
  irreg_rast <- rast(latlong[, c("xmin", "xmax", "ymin", "ymax")], 
                     values = latlong$value, 
                     crs = "EPSG:4326")
  

  该对象支持extract()等空间操作，能直接处理与多边形的重叠面积计算。

2. 转换为多要素多边形（sf包） #

将每个单元格转为多边形是最灵活的方案，适合复杂的空间叠加分析：

• 生成多边形要素：

  library(sf)
  
  # 为每行生成单元格多边形
  grid_polys <- lapply(1:nrow(latlong), function(i) {
    st_polygon(list(rbind(
      c(latlong$xmin[i], latlong$ymin[i]),
      c(latlong$xmax[i], latlong$ymin[i]),
      c(latlong$xmax[i], latlong$ymax[i]),
      c(latlong$xmin[i], latlong$ymax[i]),
      c(latlong$xmin[i], latlong$ymin[i])
    )))
  })
  
  # 转为sf对象
  grid_sf <- st_sf(value = latlong$value,
                   geometry = st_sfc(grid_polys),
                   crs = "EPSG:4326")
  

• 重叠面积计算示例：

  # 假设target_poly是你的目标多边形sf对象
  overlap <- st_intersection(grid_sf, target_poly)
  # 计算重叠面积占单元格总面积的比例
  overlap$area_ratio <- st_area(overlap) / st_area(grid_sf[overlap$row_id, ])
  # 计算加权值
  weighted_sum <- sum(overlap$value * overlap$area_ratio)
  

3. 使用stars包处理不规则网格 #

stars包原生支持不规则网格，适合多维空间数据操作：

library(stars)

# 创建stars不规则网格对象
irreg_stars <- st_as_stars(
  data.frame(
    x = latlong$x, y = latlong$y,
    xmin = latlong$xmin, xmax = latlong$xmax,
    ymin = latlong$ymin, ymax = latlong$ymax,
    value = latlong$value
  ),
  dims = st_dimensions(
    x = st_bounds(c(latlong$xmin, latlong$xmax)),
    y = st_bounds(c(latlong$ymin, latlong$ymax)),
    values = "value"
  ),
  crs = "EPSG:4326"
)

后续可通过st_extract()等函数完成与多边形的空间交互。

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