使用terra包转换XYZ数据为栅格时遇单Y坐标错误的解决求助

我用如下多边形裁剪了栅格：

我希望根据栅格单元在多边形内的面积占比调整2个单元的数值——比如单元1有90%面积在多边形内，就把它的值乘以0.9，另一单元同理。

我采用terra::extract函数，设置xy = T和weights = T提取单元值及权重，代码如下：

library(terra)
library(dplyr)
  
crop_rast_df <- terra::extract(crop_rast, my_shp , weights = T, touches = T, xy = T)
crop_rast_df
# 输出结果：
# ID GDP_PPP_3       x      y   weight
# 1  79067856 -85.875 32.125 0.874914
# 1 344945408 -85.625 32.125 0.932744
  
# 调整数值
crop_rast_df <- crop_rast_df %>% 
                  dplyr::mutate(eff_gdp = GDP_PPP_3 * weight) %>% # 考虑单元贡献占比
                  dplyr::select(x, y, eff_gdp)
crop_rast_df
# 输出结果：
# x      y   eff_gdp
# -85.875 32.125  38269172
# -85.625 32.125 177990394

但尝试将数据转回栅格时，出现如下错误：

crop_rast_adj <- terra::rast(crop_rast_df, type = 'xyz', crs = crs(temp_gdp))
# 错误信息：
# Error: [rast] cannot create a raster geometry from a single y coordinate

问题原因 #

报错的核心原因是你的crop_rast_df中只有同一个y坐标的栅格点，terra::rast(type='xyz')需要至少两组不同的x、y值来自动推断栅格的分辨率、范围等几何参数，单个y值无法构建完整的栅格结构。

修正方法 #

方法1：基于原始栅格模板更新值（推荐） #

利用原始栅格的几何参数作为模板，只替换需要调整的单元值，这是最稳妥的方式：

# 1. 复制原始栅格作为模板，保留所有单元的几何信息
crop_rast_adj <- crop_rast

# 2. 将调整后的数据转为矢量点格式，方便匹配栅格单元
adjust_points <- vect(crop_rast_df, geom = c("x", "y"), crs = crs(crop_rast))

# 3. 替换模板栅格中对应位置的数值
crop_rast_adj <- replace_values(crop_rast_adj, adjust_points, field = "eff_gdp")

方法2：手动指定栅格参数创建 #

如果一定要通过xyz数据框创建栅格，需要手动指定分辨率、范围等参数，不让函数自动推断：

# 获取原始栅格的几何参数
rast_res <- res(crop_rast)  # 分辨率
rast_ext <- ext(crop_rast)  # 范围

# 1. 创建空栅格，匹配原始栅格的几何属性
crop_rast_adj <- rast(ext = rast_ext, res = rast_res, crs = crs(crop_rast))

# 2. 将调整后的数据转为矢量点
adjust_points <- vect(crop_rast_df, geom = c("x", "y"), crs = crs(crop_rast))

# 3. 将点数据的值填充到空栅格中
crop_rast_adj <- rasterize(adjust_points, crop_rast_adj, field = "eff_gdp")

方法3：直接通过权重栅格计算（更高效） #

可以跳过extract步骤，直接生成权重栅格后与原栅格相乘，一步完成调整：

# 生成每个栅格单元在多边形内的面积占比权重栅格
weight_rast <- rasterize(my_shp, crop_rast, field = 1, touches = TRUE, background = 0)
# 注：如果需要更精确的面积占比，可使用`exactextractr`包的`exact_extract`函数

# 直接计算调整后的栅格值
crop_rast_adj <- crop_rast * weight_rast

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