MODIS8天数据转日数据：日期重复与缺失问题求解

我编写了一段代码，将MODIS 8天尺度的ET数据转换为日尺度数据。由于每年最后一期MODIS数据根据是否为闰年，需分为5天或6天而非8天，我已对此进行处理，但结果出现日期重复和缺失的问题：例如2016年9月20日的数据重复出现两次，却缺失9月21日的数据；同时时间序列的起始日期2014-01-01也不存在。

为便于理解代码逻辑，我先逐个修改了栅格的CRS（直接对栅格堆栈处理时会报错，该步骤已成功完成），随后基于重投影后的栅格列表进行循环处理，问题出在该循环部分。

以下是我的代码：

library(raster)
library(terra)
library(lubridate)

archivos_ET <- list.files(path = "/", pattern = "ET_\\d{4}_\\d{3}", full.names = TRUE)

#Save all raster et files in a list called et_proy

# Create a list to store the daily rasters
et_dia <- list()

for (i in 1:length(et_proy)) {
  # # Load the raster
  raster_actual <- et_proy[[i]]
  # Get the year and day of the current file name
  partes_nombre <- strsplit(basename(archivos_ET[i]), "_")[[1]]
  ano <- as.integer(partes_nombre[2])
  dia <- as.integer(gsub(".tif", "", partes_nombre[3]))
  # Check if the year is a leap year
  es_bisiesto <- ifelse(ano %% 4 == 0 & (ano %% 100 != 0 | ano %% 400 == 0), TRUE, FALSE)
  # Determine the divisor according to whether it is a leap year and the day is 361 (it is the last day of the modis products in a year)
  divisor <- ifelse(dia == 361, ifelse(es_bisiesto, 6, 5), 8)
  # Loop to split the current raster into daily rasters and store them
  for (j in 1:(divisor-1)) {
    # Calcular la fecha del raster diario
    fecha <- as.Date(paste(ano, "-01-01", sep = "")) + dia + j - 1
    # Calculate date from daily raster
    raster_diario <- raster_actual / divisor
    # Assign date to daily raster
    names(raster_diario) <- format(fecha, "%Y-%m-%d")
    # Add the daily raster to the list
    et_dia[[length(et_dia) + 1]] <- raster_diario
  }
  # Calculate the date of the last raster
  fecha_ultimo <- as.Date(paste(ano, "-01-01", sep = "")) + dia + divisor - 2
  # Calculate the last raster
  raster_ultimo <- raster_actual / divisor
  # Assign date to last raster
  names(raster_ultimo) <- format(fecha_ultimo, "%Y-%m-%d")
  # Add the last raster to the list
  et_dia[[length(et_dia) + 1]] <- raster_ultimo
}

结果示例如下：

[[993]]
class      : RasterLayer 
dimensions : 33, 56, 1848  (nrow, ncol, ncell)
resolution : 346, 463  (x, y)
extent     : 602219.9, 621595.9, 5353158, 5368437  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs 
source     : memory
names      : X2016.09.20 
values     : 0.6063787, 2.016294  (min, max)


[[994]]
class      : RasterLayer 
dimensions : 33, 56, 1848  (nrow, ncol, ncell)
resolution : 346, 463  (x, y)
extent     : 602219.9, 621595.9, 5353158, 5368437  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs 
source     : memory
names      : X2016.09.20 
values     : 0.6063787, 2.016294  (min, max)


[[995]]
class      : RasterLayer 
dimensions : 33, 56, 1848  (nrow, ncol, ncell)
resolution : 346, 463  (x, y)
extent     : 602219.9, 621595.9, 5353158, 5368437  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs 
source     : memory
names      : X2016.09.22 
values     : 1.3125, 2.177313  (min, max)


[[996]]
class      : RasterLayer 
dimensions : 33, 56, 1848  (nrow, ncol, ncell)
resolution : 346, 463  (x, y)
extent     : 602219.9, 621595.9, 5353158, 5368437  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs 
source     : memory
names      : X2016.09.23 
values     : 1.3125, 2.177313  (min, max)

此外，我还发现时间序列的起始数据2014-01-01也缺失了。

问题根源分析 #

代码的核心错误出在日期计算逻辑：

1. 日期偏移计算错误：MODIS产品的儒略日是周期起始日，但代码中用dia + j -1计算日期，直接跳过了起始日，且后续偏移导致日期错位。
2. 重复添加日期：内层循环和单独添加的最后一个日期计算重复，导致同一日期被存入两次，同时跳过了中间的日期。

修正后的代码 #

library(raster)
library(terra)
library(lubridate)

archivos_ET <- list.files(path = "/", pattern = "ET_\\d{4}_\\d{3}", full.names = TRUE)

# 假设et_proy已经是重投影后的栅格列表
et_dia <- list()

for (i in 1:length(et_proy)) {
  raster_actual <- et_proy[[i]]
  # 解析文件名中的年份和儒略日
  partes_nombre <- strsplit(basename(archivos_ET[i]), "_")[[1]]
  ano <- as.integer(partes_nombre[2])
  dia_juliano <- as.integer(gsub(".tif", "", partes_nombre[3]))
  
  # 用lubridate内置函数判断闰年，更准确可靠
  es_bisiesto <- leap_year(ano)
  
  # 确定当前周期的天数
  divisor <- if (dia_juliano == 361) {
    ifelse(es_bisiesto, 6, 5)
  } else {
    8
  }
  
  # 生成当前周期的连续日期序列
  fecha_inicio <- as.Date(dia_juliano, origin = paste0(ano, "-01-01"))
  fechas <- seq(fecha_inicio, by = "day", length.out = divisor)
  
  # 拆分栅格并添加到列表
  for (fecha in fechas) {
    raster_diario <- raster_actual / divisor
    names(raster_diario) <- format(fecha, "%Y-%m-%d")
    et_dia[[length(et_dia) + 1]] <- raster_diario
  }
}

修正说明 #

1. 简化闰年判断：使用lubridate::leap_year替代手动计算，避免逻辑错误。
2. 正确生成日期序列：从MODIS产品的起始儒略日开始，直接生成连续的divisor个日期，确保起始日期不丢失、后续日期无错位。
3. 合并循环逻辑：用单个循环遍历日期序列，避免重复添加和日期缺失问题。

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