How to make the result of dates in the loop correlative?

Question

I made a code to calculate daily data from MODIS 8-day data, as it is known the last data of the year is always divided by 5 or 6 and not by 8 depending on whether it is a leap year or not so I added it, the problem is the result because I have repeated dates and other missing ones, for example September 20 is repeated twice but I don't have 21. How could I solve this?

So that the code is well understood, the first thing I did was change the crs of my rasters one by one because I had errors applying it directly to the stack at the end, this worked fine, and then with that list with my reprojected rasters I worked with the loop I have doubts

Here is my code

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
}

And here is an example of my result

[[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)

I also realized that the first data of my series of years, that is, 2014-01-01, is not there either.

I’m voting to close. This should be migrated to the SE Geographic Information Systems forum — IRTFM, May 14 '23 at 20:54
I’m puzzled. I thought MODIS data was satellite imagery. You’re working with data that is encapsulated with coordinate reference systems. That meant to me that the GIS forum would be the place to find people who deal with that on a regular basis. What am I reading wrong? If I’m wrong then describe the problem more fully and include both data and `library` calls to load all needed packages. — IRTFM, May 14 '23 at 22:25
@IRTFM I'm sorry, maybe my question was not well understood but it doesn't have to do so much with MODIS, rather with the logic of the loop, I don't know what I'm doing wrong because the result in relation to their names (dates) is not as expected — Vanesa Palma, May 15 '23 at 03:42
As I said. No data and incomplete code …. No possibility of testing or solid understanding. — IRTFM, May 15 '23 at 14:49
I think the question if suitable in principle, but can you please improve it. First, make it *minimal* (only include packages needed (never "library(tidyverse)"), and no pre-processing such as projectRaster. Also make it *self contained* by adding example data (create a raster with code that has the same number of layers / dates as your data; but only a few cells to keep things light). In fact, you may not need **any** raster data. To start with you could show what you want with a simple vector. Show input and desired output. Once that is solved, the raster bit is probably trivial. — Robert Hijmans, May 15 '23 at 16:33

Robert Hijmans · Answer 1 · 2023-05-15T20:20:14.507

Here is a general solution.

Example data at 8-day time step

library(terra)
r <- rast(system.file("ex/logo.tif", package="terra"))   
s <- c(r/c(1:3), r/c(4:6))
time(s) <- as.Date("2001-01-03") + seq(0, 40, 8)
time(s)
#[1] "2001-01-03" "2001-01-11" "2001-01-19" "2001-01-27" "2001-02-04" "2001-02-12"

Add one layer to make the smallest time-step one day

x <- s[[1]]
time(x) <- time(x)-1
x <- c(x, s)

Now use fillTime to add (empty) layers inbetween the 8 days intervals, and remove the added layer.

sx <- fillTime(x)[[-1]]
time(sx)
# [1] "2001-01-03" "2001-01-04" "2001-01-05" "2001-01-06" "2001-01-07" "2001-01-08" "2001-01-09" "2001-01-10" "2001-01-11" "2001-01-12" "2001-01-13"
#[12] "2001-01-14" "2001-01-15" "2001-01-16" "2001-01-17" "2001-01-18" "2001-01-19" "2001-01-20" "2001-01-21" "2001-01-22" "2001-01-23" "2001-01-24"
#[23] "2001-01-25" "2001-01-26" "2001-01-27" "2001-01-28" "2001-01-29" "2001-01-30" "2001-01-31" "2001-02-01" "2001-02-02" "2001-02-03" "2001-02-04"
#[34] "2001-02-05" "2001-02-06" "2001-02-07" "2001-02-08" "2001-02-09" "2001-02-10" "2001-02-11" "2001-02-12" "2001-02-13" "2001-02-14" "2001-02-15"
#[45] "2001-02-16" "2001-02-17" "2001-02-18" "2001-02-19" "2001-02-20" "2001-02-21" "2001-02-22" "2001-02-23" "2001-02-24" "2001-02-25" "2001-02-26"
#[56] "2001-02-27" "2001-02-28" "2001-03-01" "2001-03-02" "2001-03-03" "2001-03-04" "2001-03-05" "2001-03-06" "2001-03-07" "2001-03-08"

Use approximate to estimate missing values by linear interpolation.

a <- approximate(sx)

plot(time(a), minmax(a)[2,])

A more "manual" approach:

stm <- time(s)
tms <- seq(stm[1], stm[length(stm)], 1)
r <- rast(s, nlyr=length(tms))
time(r) <- tms

i <- match(stm, tms)
r[[i]] <- s

And take it from there.

Thank you very much for answering, I will try it with this method! — Vanesa Palma, May 16 '23 at 01:53

score 0 · Answer 2 · answered May 16 '23 at 02:08

I think I have solved my question, the solution was to have created an initial date and to have reviewed where the name was assigned.


et_dia <- list()
fecha_actual <- as.Date("2013-12-31")

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 <- as.integer(substr(partes_nombre[3], 1, 3))
  es_bisiesto <- ifelse(ano %% 4 == 0 & (ano %% 100 != 0 | ano %% 400 == 0), TRUE, FALSE)
  divisor <- ifelse(dia == 361, ifelse(es_bisiesto, 6, 5), 8)
  for (j in 1:(divisor)) {
    fecha_actual <- fecha_actual + 1
    raster_diario <- raster_actual / divisor
    names(raster_diario) <- format(fecha_actual, "%Y-%m-%d")
    et_dia[[length(et_dia) + 1]] <- raster_diario
  }
}

How to make the result of dates in the loop correlative?

2 Answers2