R语言实现Sentinel-2影像插值：缺失值修复效果不稳定求助

我有连续数月的Sentinel-2影像数据，尝试用邻近日期的影像值对缺失值进行插值填补，但当前的R代码效果不稳定——部分缺失值能被成功替换，另一部分则无法填充。

原尝试代码 #

library(raster)
library(foreach)
library(doParallel)

cl <- makeCluster(30)
registerDoParallel(cl)

path <- "/R//umixing//havel_sentinel2//"
filename <- list.files(path)

df_list <- foreach(rfile = filename, .combine = cbind, .packages = c("raster")) %dopar% {
  r <- stack(paste0(path, rfile))
  ras <- round(r, 4)
  df1 <- as.data.frame(ras, xy=TRUE)
  drops <- c("x","y")
  df1 <- df1[ , !(names(df1) %in% drops)]
  colnames(df1) <- c("band1", "band2","band3", "band4", "band5", "band6", "band7", "band8", "band9", "band10")
  
  return (df1)
}

x_parallel <- function(raster, band_count = 10) {
  
  sentinel <- -3.2768
  tol      <- 1e-6
  
  # Helper to find the next non‐NA in a given direction
  findNearestNonNA <- function(row_data, col_idx, step) {
    new_col <- col_idx + step
    while (
      new_col > 0 && new_col <= length(row_data) &&
      (is.na(row_data[new_col]) ||
       abs(row_data[new_col] - sentinel) < tol)
    ) {
      new_col <- new_col + step
    }
    if (
      new_col <= 0 ||
      new_col > length(row_data) ||
      is.na(row_data[new_col]) ||
      abs(row_data[new_col] - sentinel) < tol
    ) {
      return(NA)
    } else {
      return(row_data[new_col])
    }
  }
  
  # Set up parallel back end
  cores <- detectCores()
  cl <- makeCluster(30) 
  registerDoParallel(cl)
  
  # Identify only the rows that need fixing
  rows_to_fix <- which(apply(raster, 1, function(row) any(is.na(row) | abs(row - (-3.2768)) < 1e-6)))
  
  raster_fixed  <- raster
  # Process each row in parallel
  fixed_rows <- foreach(r = rows_to_fix, .combine=rbind) %dopar% {
    row_data <- raster[r, ]
    na_cols <- which(is.na(row_data) | abs(row_data - (-3.2768)) < 1e-6)
    if (length(na_cols) == 0) {
      # Nothing to fix
      row_data
    } else {
      for (c_na in na_cols) {
        val1 <- findNearestNonNA(row_data, c_na, +band_count)
        val2 <- findNearestNonNA(row_data, c_na, -band_count)
        if (!is.na(val1) && !is.na(val2)) {
          row_data[c_na] <- mean(c(val1, val2))
        } else if (!is.na(val1)) {
          row_data[c_na] <- val1
        } else if (!is.na(val2)) {
          row_data[c_na] <- val2
        }
      }
      row_data
    }
  }
  
  stopCluster(cl)
  # Put the fixed rows back into the full raster
  raster_fixed[rows_to_fix, ] <- fixed_rows # Clean up
  return(raster_fixed)
}

# 1) Run the interpolation on the combined dataframe
df_interpolated <- x_parallel(df_list, band_count = 10)


# 2) Split into single images (10 columns per image) and save
for(i in seq_along(filename)) {
  
  # Read original file so we can match its raster shape & coordinate info
  base_ras <- stack(file.path(path, filename[i]))
  
  # We'll re-extract x,y from 'base_ras' so we can rebuild a RasterBrick
  df_xy    <- as.data.frame(base_ras, xy = TRUE)[, c("x","y")]
  
  # Select the 10 columns corresponding to this image
  start_col <- (i-1)*10 + 1
  end_col   <- i*10
  sub_df    <- df_interpolated[, start_col:end_col]
  
  # Name the bands so the output has clear layer names
  colnames(sub_df) <- paste0("band", 1:10)
  
  # Combine x, y with our 10-band data
  final_df <- cbind(df_xy, sub_df)
  
  # Convert to a RasterBrick with 10 layers
  rbrick <- rasterFromXYZ(final_df, crs = raster::crs(base_ras))
  
  # Write out a new file (change name/path as desired)
  out_name <- paste0("interpolated_", filename[i])
  writeRaster(rbrick, filename = file.path(path, out_name), overwrite=TRUE)
}

核心问题排查 #

1. 缺失值判断逻辑漏洞
   原代码同时判断NA和接近-3.2768的值，但浮点数精度问题可能导致部分缺失值未被识别；若数据中存在其他无效值（如未缩放的原始-32768），也会被遗漏。

2. 插值逻辑局限性
   当前仅查找缺失值前后最近的单个非缺失值，若该方向所有同波段日期值均缺失，则无法填充；且未考虑日期时间间隔，插值结果合理性不足。

3. 数据结构与内存问题
   将所有影像波段绑定为超大DataFrame，影像数量较多时会导致内存溢出或计算中断，表现为部分缺失值未处理。

4. 并行处理冲突
   外层已创建并行集群，x_parallel函数内部再次创建新集群，易引发资源冲突，导致部分任务执行失败。

修复方案与优化代码 #

优化思路 #

• 统一缺失值处理：将-3.2768转为NA，简化判断逻辑。
• 时间序列线性插值：基于日期维度对每个像素-波段进行插值，利用时间连续性提升效果。
• 内存优化：按波段处理时间序列，避免超大DataFrame。
• 简化并行逻辑：避免嵌套集群，改用按波段并行。

修改后的代码示例 #

library(raster)
library(tidyverse)
library(furrr)
library(lubridate)

# 设置并行
plan(multisession, workers = 10)

path <- "/R//umixing//havel_sentinel2//"
filename <- list.files(path, full.names = TRUE)

# 从文件名提取日期（假设文件名含8位日期，如"S2A_20230101.tif"）
extract_date <- function(file) {
  str_extract(basename(file), "\\d{8}") %>% ymd()
}

dates <- map_dbl(filename, extract_date) %>% sort()
sorted_files <- filename[order(dates)]

# 读取所有影像为时间序列栈
ts_stack <- stack(sorted_files)
names(ts_stack) <- paste0("date_", str_remove_all(dates, "-"))

# 统一缺失值为NA
ts_stack[ts_stack == -3.2768] <- NA

# 单波段时间插值函数
interpolate_band <- function(band_raster) {
  ts_df <- as.data.frame(band_raster, xy = TRUE) %>%
    pivot_longer(cols = starts_with("date_"), names_to = "date", values_to = "value") %>%
    mutate(date = str_remove(date, "date_") %>% ymd()) %>%
    group_by(x, y) %>%
    # 线性插值，首尾缺失值用最近值填充
    mutate(value = na.approx(value, x = date, na.rm = FALSE, rule = 2)) %>%
    ungroup() %>%
    pivot_wider(names_from = "date", values_from = "value")
  
  raster_from_xyz(ts_df, crs = crs(band_raster))
}

# 并行插值所有波段
interpolated_stack <- map(1:nlayers(ts_stack), function(i) {
  interpolate_band(ts_stack[[i]])
}) %>% stack()

# 拆分保存为单个影像
walk(seq_along(sorted_files), function(i) {
  date_str <- str_remove_all(dates[i], "-")
  output_file <- file.path(path, paste0("interpolated_S2_", date_str, ".tif"))
  writeRaster(interpolated_stack[[i]], filename = output_file, overwrite = TRUE)
})

# 关闭并行
plan(sequential)

关键改进说明 #

1. 统一缺失值：直接替换-3.2768为NA，避免精度判断误差。
2. 时间线性插值：使用na.approx结合日期维度插值，rule=2处理首尾缺失值。
3. 内存优化：按波段处理时间序列，降低内存占用。
4. 并行简化：用furrr并行映射，避免嵌套集群冲突。

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