如何生成满足约束且可参考分布的5个随机raster？

一、生成满足数值范围与总和为1的随机栅格 #

以下针对两种常见的“总和为1”需求分别提供解决方案：

情况1：每个栅格自身的所有像素值之和为1 #

每个栅格的单个像素值在指定xmin[i]到xmax[i]范围内，且整个栅格的像素总和为1。使用R的raster包结合约束优化实现：

library(raster)

# 定义栅格尺寸（可根据需求修改）
nrow <- 10
ncol <- 10

# 用户给定的数值范围
xmin <- c(0, 0, 0, 0, 0)
xmax <- c(26.0, 78.6, 14.4, 39.4, 70.8)

# 生成单个约束栅格的函数
generate_constrained_raster <- function(min_val, max_val, nrow, ncol, target_sum = 1) {
  ncells <- nrow * ncol
  # 生成初始随机值
  init_vals <- runif(ncells, min = min_val, max = max_val)
  
  # 若初始总和已符合要求，直接返回
  if (abs(sum(init_vals) - target_sum) < 1e-6) {
    return(raster(nrow = nrow, ncol = ncol, vals = init_vals))
  }
  
  # 约束优化：调整值使总和为target_sum，同时保持在范围内
  objective <- function(vals) abs(sum(vals) - target_sum)
  result <- optim(init_vals, objective, method = "L-BFGS-B",
                  lower = rep(min_val, ncells), upper = rep(max_val, ncells),
                  control = list(maxit = 1000))
  
  raster(nrow = nrow, ncol = ncol, vals = result$par)
}

# 批量生成5个栅格
single_sum_rasters <- lapply(1:5, function(i) {
  generate_constrained_raster(xmin[i], xmax[i], nrow, ncol)
})

# 验证结果
lapply(single_sum_rasters, function(r) {
  cat("栅格总和：", round(sum(values(r)), 6), "\n")
  cat("像素范围：", round(min(values(r)), 2), "~", round(max(values(r)), 2), "\n\n")
})

情况2：同一位置的5个栅格像素值之和为1 #

每个对应位置的5个栅格像素值相加等于1，且每个值在各自的xmin[i]到xmax[i]范围内：

# 生成联合约束栅格的函数
generate_joint_rasters <- function(xmin, xmax, nrow, ncol) {
  ncells <- nrow * ncol
  n_rasters <- length(xmin)
  vals_matrix <- matrix(0, nrow = ncells, ncol = n_rasters)
  
  # 对每个像素位置生成符合条件的5个值
  for (i in 1:ncells) {
    init_vals <- runif(n_rasters, min = xmin, max = xmax)
    objective <- function(vals) abs(sum(vals) - 1)
    result <- optim(init_vals, objective, method = "L-BFGS-B",
                    lower = xmin, upper = xmax,
                    control = list(maxit = 1000))
    vals_matrix[i, ] <- result$par
  }
  
  # 转换为栅格列表
  lapply(1:n_rasters, function(i) {
    raster(nrow = nrow, ncol = ncol, vals = vals_matrix[, i])
  })
}

# 生成5个联合约束栅格
joint_sum_rasters <- generate_joint_rasters(xmin, xmax, nrow, ncol)

# 验证随机位置的和
sample_cells <- sample(ncell(joint_sum_rasters[[1]]), 5)
for (cell in sample_cells) {
  cell_sum <- sum(sapply(joint_sum_rasters, function(r) values(r)[cell]))
  cat("像素位置", cell, "的和：", round(cell_sum, 6), "\n")
}

二、让随机栅格匹配已有栅格的数值分布 #

使用分位数匹配方法，让生成的栅格与参考栅格拥有相同的分布特征，同时满足范围和总和约束：

# 示例参考栅格（替换为你的实际参考栅格）
ref_raster <- raster(nrow = nrow, ncol = ncol)
values(ref_raster) <- rnorm(ncell(ref_raster), mean = 40, sd = 15) # 模拟参考分布

# 生成匹配分布的约束栅格函数
generate_matched_raster <- function(min_val, max_val, nrow, ncol, ref_raster, target_sum = 1) {
  ncells <- nrow * ncol
  # 生成初始范围内随机值
  init_vals <- runif(ncells, min = min_val, max = max_val)
  
  # 分位数匹配：将初始值的分布映射到参考栅格的分布
  ref_quant <- quantile(values(ref_raster), probs = seq(0, 1, length.out = ncells))
  init_quant <- quantile(init_vals, probs = seq(0, 1, length.out = ncells))
  matched_vals <- approx(init_quant, ref_quant, xout = init_vals)$y
  
  # 截断到指定范围，缩放总和到1
  matched_vals <- pmax(min_val, pmin(max_val, matched_vals))
  scale_factor <- target_sum / sum(matched_vals)
  final_vals <- pmax(min_val, pmin(max_val, matched_vals * scale_factor))
  
  raster(nrow = nrow, ncol = ncol, vals = final_vals)
}

# 批量生成5个匹配分布的栅格
matched_rasters <- lapply(1:5, function(i) {
  generate_matched_raster(xmin[i], xmax[i], nrow, ncol, ref_raster)
})

# 验证分布匹配（用QQ图可视化）
qqplot(values(ref_raster), values(matched_rasters[[1]]), 
       xlab = "参考栅格值", ylab = "生成栅格值", main = "分位数匹配验证")
abline(0, 1, col = "red", lwd = 2)

注意事项 #

• 如果参考栅格的分布范围与目标范围差异过大，分位数匹配后可能需要截断，这会轻微改变分布特征，建议提前确保参考分布的核心范围在目标范围内。
• 优化算法的maxit参数可根据需求调整，确保收敛到符合要求的结果。

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