考虑NA像素的多栅格平均：MODIS LST数据计算异常排查

我使用MODIS11A2的8天平均LST栅格（.tif）计算年度日均温，2001年共44个文件（缺失2个）。需求如下：

• 仅当某一像素所有栅格均为NA时，才排除该像素的平均计算
• 若像素仅部分栅格为NA（比如1个NA、43个有效数据），则用有效数据总和除以有效数量计算平均值

但当前代码会把这类仅部分NA的像素也排除，求排查解决。

原代码 #

# Load necessary libraries
library(raster)
library(rgdal)

# Set the path to your folder containing the .tif files
folder_path <- "D:/A/LST/mosaic_clipped_reprojected/test"

# List all .tif files in the folder
tif_files <- list.files(folder_path, pattern = "*.tif$", full.names = TRUE)

# Initialize variables to store the sum of rasters and valid count
total_sum <- NULL
valid_count <- NULL

# Set NoData value for MODIS LST data
nodata_value <- -9999

# Loop through all the tif files and process them
for (tif_file in tif_files) {
  # Read the raster file
  raster_data <- raster(tif_file)

  # Set the NoData value for the raster (this step ensures NoData values are handled correctly)
  raster_data[raster_data == nodata_value] <- NA  # Convert NoData to NA

  # Initialize total_sum and valid_count on the first iteration
  if (is.null(total_sum)) {
    total_sum <- raster_data
    valid_count <- !is.na(raster_data)  # Create a mask for valid data (TRUE = valid, FALSE = NA)
  } else {
    # Add current raster data to the total_sum, and count valid pixels
    total_sum <- total_sum + raster_data
    valid_count <- valid_count + !is.na(raster_data)  # Count valid pixels (add 1 for valid, 0 for NA)
  }
}

# After all rasters are processed, replace invalid count (0 values) with NA to avoid dividing by 0
valid_count[valid_count == 0] <- NA

# Calculate the daily average by dividing the total sum by the valid count for each pixel
# This ensures that only valid pixels contribute to the sum and the division.
daily_average <- total_sum / valid_count

# Define the reference raster for resampling (choose one of your original rasters)
reference_raster <- raster(tif_files[1])  # Example: using the first raster file as reference

# Resample the daily average to the reference raster using bilinear interpolation
resampled_daily_average <- resample(daily_average, reference_raster, method = "bilinear")

# Define the output file path for the daily average raster
output_file <- "D:/A/LST/AnnualAvgLST/daily_average_resampled_test1.tif"

# Write the resampled daily average raster to a new file
writeRaster(resampled_daily_average, output_file, format = "GTiff", overwrite = TRUE)

# Print message to indicate the result
cat("Resampled daily average raster has been saved to:", output_file, "\n")

问题根源 #

核心问题是栅格求和时的NA传播：R中NA + 数值 = NA，如果某像素在第一个栅格中是NA，后续即使有有效数据，累加后的总和仍会保持NA，最终计算平均值时该像素会被判定为无效，不符合需求。

修正后的代码 #

# Load necessary libraries
library(raster)
library(rgdal)

# Set the path to your folder containing the .tif files
folder_path <- "D:/A/LST/mosaic_clipped_reprojected/test"

# List all .tif files in the folder
tif_files <- list.files(folder_path, pattern = "*.tif$", full.names = TRUE)

# Initialize variables to store the sum of rasters and valid count
total_sum <- NULL
valid_count <- NULL

# Set NoData value for MODIS LST data
nodata_value <- -9999

# Loop through all the tif files and process them
for (tif_file in tif_files) {
  # Read the raster file
  raster_data <- raster(tif_file)
  
  # Convert NoData to NA, then create a version where NA becomes 0 for summation
  raster_data[raster_data == nodata_value] <- NA
  raster_sum <- raster_data
  raster_sum[is.na(raster_sum)] <- 0  # 替换NA为0，避免求和时传播NA
  
  # Initialize total_sum and valid_count on the first iteration
  if (is.null(total_sum)) {
    total_sum <- raster_sum
    valid_count <- !is.na(raster_data)  # 基于原始数据统计有效像素
  } else {
    # 累加有效数据的数值，累加有效计数
    total_sum <- total_sum + raster_sum
    valid_count <- valid_count + !is.na(raster_data)
  }
}

# 将完全无数据的像素设为NA
valid_count[valid_count == 0] <- NA

# 计算平均值：总和除以有效计数
daily_average <- total_sum / valid_count

# 重采样（若原始栅格分辨率/范围一致，可省略此步骤）
reference_raster <- raster(tif_files[1])
resampled_daily_average <- resample(daily_average, reference_raster, method = "bilinear")

# 输出结果
output_file <- "D:/A/LST/AnnualAvgLST/daily_average_resampled_test1.tif"
writeRaster(resampled_daily_average, output_file, format = "GTiff", overwrite = TRUE)

cat("Resampled daily average raster has been saved to:", output_file, "\n")

关键修改说明 #

• 创建raster_sum临时对象，将NA替换为0后再参与求和，彻底避免NA在累加过程中传播
• valid_count仍基于原始raster_data的NA状态统计，确保有效数据的计数准确
• 最终仅将有效计数为0（完全无数据）的像素设为NA，符合需求

内容的提问来源于stack exchange，提问作者Alexia k Boston