英文:
How to make the result of dates in the loop correlative?
问题
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 rasterset_dia <- list()for (i in 1:length(et_proy)) {# Load the rasterraster_actual <- et_proy[[i]]# Get the year and day of the current file namepartes_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 yeares_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 themfor (j in 1:(divisor-1)) {# Calculate the date of the daily rasterfecha <- as.Date(paste(ano, "-01-01", sep = "")) + dia + j - 1# Calculate the daily rasterraster_diario <- raster_actual / divisor# Assign date to daily rasternames(raster_diario) <- format(fecha, "%Y-%m-%d")# Add the daily raster to the listet_dia[[length(et_dia) + 1]] <- raster_diario}# Calculate the date of the last rasterfecha_ultimo <- as.Date(paste(ano, "-01-01", sep = "")) + dia + divisor - 2# Calculate the last rasterraster_ultimo <- raster_actual / divisor# Assign date to last rasternames(raster_ultimo) <- format(fecha_ultimo, "%Y-%m-%d")# Add the last raster to the listet_dia[[length(et_dia) + 1]] <- raster_ultimo}
And here is an example of my result
[[993]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.20values : 0.6063787, 2.016294 (min, max)[[994]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.20values : 0.6063787, 2.016294 (min, max)[[995]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.22values : 1.3125, 2.177313 (min, max)[[996]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.23values : 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 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 rasterset_dia <- list()for (i in 1:length(et_proy)) {# # Load the rasterraster_actual <- et_proy[[i]]# Get the year and day of the current file namepartes_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 yeares_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 themfor (j in 1:(divisor-1)) {# Calcular la fecha del raster diariofecha <- as.Date(paste(ano, "-01-01", sep = "")) + dia + j - 1# Calculate date from daily rasterraster_diario <- raster_actual / divisor# Assign date to daily rasternames(raster_diario) <- format(fecha, "%Y-%m-%d")# Add the daily raster to the listet_dia[[length(et_dia) + 1]] <- raster_diario}# Calculate the date of the last rasterfecha_ultimo <- as.Date(paste(ano, "-01-01", sep = "")) + dia + divisor - 2# Calculate the last rasterraster_ultimo <- raster_actual / divisor# Assign date to last rasternames(raster_ultimo) <- format(fecha_ultimo, "%Y-%m-%d")# Add the last raster to the listet_dia[[length(et_dia) + 1]] <- raster_ultimo}
And here is an example of my result
[[993]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.20values : 0.6063787, 2.016294 (min, max)[[994]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.20values : 0.6063787, 2.016294 (min, max)[[995]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.22values : 1.3125, 2.177313 (min, max)[[996]]class : RasterLayerdimensions : 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_defssource : memorynames : X2016.09.23values : 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.
答案1
得分: 1
以下是翻译好的部分：
这是一个通用解决方案。
8天时间步骤的示例数据
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"
添加一个层以使最小时间步长为一天
x <- s[[1]]time(x) <- time(x)-1x <- c(x, s)
现在使用fillTime在8天间隔之间添加（空）层，并删除添加的层。
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"
使用approximate进行线性插值估算缺失值。
a <- approximate(sx)
plot(time(a), minmax(a)[2,])
更“手动”的方法：
stm <- time(s)tms <- seq(stm[1], stm[length(stm)], 1)r <- rast(s, nlyr=length(tms))time(r) <- tmsi <- match(stm, tms)r[[i]] <- s
然后继续处理。
英文:
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)-1x <- 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) <- tmsi <- match(stm, tms)r[[i]] <- s
And take it from there.
答案2
得分: 0
我认为我已经解决了我的问题，解决方法是创建了一个初始日期并检查了名称的分配位置。
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 + 1raster_diario <- raster_actual / divisornames(raster_diario) <- format(fecha_actual, "%Y-%m-%d")et_dia[[length(et_dia) + 1]] <- raster_diario}}
英文:
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 + 1raster_diario <- raster_actual / divisornames(raster_diario) <- format(fecha_actual, "%Y-%m-%d")et_dia[[length(et_dia) + 1]] <- raster_diario}}
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论