Group raster files by week/month

Question

I have downloaded AOD data ("MCD19A2") using the MODIStsp package for several years.

library(raster)
library(stars)

The data I have extracted has the following format:

MCD19A2_Optical_Depth_055_2011_001.tif
MCD19A2_Optical_Depth_055_2011_002.tif ...
MCD19A2_Optical_Depth_055_2011_365.tif

Where files from each year are in a different folder.

I would like to be able to have weekly and monthly means and vectorize them. I have done this before with all the files from a year folder this way:

path_21<-paste0(wd,"/2021/Optical_Depth_055")
#list all files from the 2021 folder
all_modis_21<-list.files(path_21,
                         full.names = TRUE,
                         pattern = "*.tif$")

all_modis_21_st <- stack(all_modis_21)
all_modis_21_br<-brick(all_modis_21_st)

# calculate the mean of all files
mean_2021<-calc(all_modis_21_br, fun = mean, na.rm = TRUE)
path_21_mean<-paste0(wd,"MEANS/2021")
# store the mean in the "MEANS" folder
writeRaster(mean_2021, filename=file.path(path_21_mean,"mean_2021.tiff"),
            format="GTiff", overwrite=TRUE)
mean_2021<-read_stars("Data/Pollution/MODIS/NASA_DATA/MEANS/2021/mean_2021.tiff")
wrap21 = st_warp(mean_2021, crs = 25831)
# vectorize the mean for 2021
p = st_as_sf(wrap21)

The previous code worked for creating the mean and vectorizing the files for the whole folder. But I don't know how to group them by week.

So far I have changed the name of the files to be able to work with dates.

to<-paste(as.Date(substr(list.files(path_21,pattern = "*.tif$"),start = 27,stop = 34),"%Y_%j"),".tif", sep="")
from<-list.files(path_21,pattern = "*.tif$")

raw_21<-data_frame(filename<-list.files(path_21,
                                        pattern = "*.tif$"))

And the structure of the data for each folder looks like this:

head(raw_21)
 
1 2011-01-01.tif                                       
 2 2011-01-02.tif                                       
 3 2011-01-03.tif                                       
 4 2011-01-04.tif                                       
 5 2011-01-05.tif                                       
 6 2011-01-06.tif                                       
 7 2011-01-07.tif                                       
 8 2011-01-08.tif                                       
 9 2011-01-09.tif                                       
10 2011-01-10.tif

But I am fairly new to R and I am lost on how to proceed from here. I would appreciate any help on how to follow :)

Please read the information at the top of the [tag:r] tag and, in particular, make your question reproducible including all inputs (using `dput`) and `library statements so anyone else can copy the code and input from the question and paste it into their session. We don't have your exact directory structure so the `setwd` should not be in the question. — G. Grothendieck, Oct 15 '22 at 13:48
It is not clear what you mean with `vectorize` in this context. You should not need the filenames, but rather extract the dates, and use that to group. But have a look at `terra::tapp`. Please do include some example data (as in the R help files), but do not use `dput`, that does not work well with spatial data. — Robert Hijmans, Oct 15 '22 at 18:20
Hi, I modified the description to make it more understandable. Does it make more sense now? — Berta_94, Oct 16 '22 at 16:08

Robert Hijmans · Accepted Answer · 2022-10-17T16:50:18.660

You want to summarize daily raster data by week, project to a different CRS, and then transform the raster data to polygons.

Here are some example filenames:

ff <- paste0("MCD19A2_Optical_Depth_055_2011_", c("001", "002", "040", "090", 120:125, 200:205, 300), ".tif")
head(ff)
#[1] "MCD19A2_Optical_Depth_055_2011_001.tif"
#[2] "MCD19A2_Optical_Depth_055_2011_002.tif"
#[3] "MCD19A2_Optical_Depth_055_2011_040.tif"
#[4] "MCD19A2_Optical_Depth_055_2011_090.tif"
#[5] "MCD19A2_Optical_Depth_055_2011_120.tif"
#[6] "MCD19A2_Optical_Depth_055_2011_121.tif"

Now you can extract the date like this

d <- gsub("MCD19A2_Optical_Depth_055_", "", basename(ff))
year <- as.numeric(substr(d, 1, 4))
doy  <- as.numeric(substr(d, 6, 9))
date <- as.Date(doy, origin=paste(year-1, "-12-31", sep=""))

Create a SpatRaster, set the dates, and get the mean for each week

x <- rast(ff)
terra::time(x) <- date
r <- tapp(x, "week", mean, na.rm=TRUE)

But note that in this year, Jan 1 is part of week 52. See ?strftime.

week <- strftime(date, format = "%V")
week
# [1] "52" "52" "06" "13" "17" "17" "18" "18" "18" "18" "29" "29"
# [13] "29" "29" "29" "29" "43"

You can also define your own 7-day periods that start on Jan 1 like this

myweek <- floor(doy / 7) + 1
myweek <- formatC(myweek, width=2, flag=0)
myweek
#[1] "01" "01" "06" "13" "18" "18" "18" "18" "18" "18" "29" "29" "29" "30" "30" "30" "43"

And use that with tapp

r <- tapp(x, myweek, mean, na.rm=TRUE)

Now you can project and create polygons

r <- project(r, "epsg:25831")
p <- as.polygons(r)
writeVector(p, "MCD19A2_Optical_Depth.shp")

A more "manual" approach could be to write a loop

uweeks <- unique(week)
rlst <- lapply(uweeks, function(w) {
    ffw <- ff[week == w]
    mean(rast(ffw))
}
r <- rast(rlst)
r <- project(r, "epsg:25831")
p <- as.polygons(r)

It works the same for months. To get the month from a date, you can do

format(date,"%m")
#[1] "01" "01" "02" "03" "04" "05" "05" "05" "05" "05" "07" "07" "07" "07" "07" "07" "10"

Final notes:

Transforming raster data to polygons is not a good idea. It is almost certainly a major mistake in your workflow. But as you are not saying why you are doing this, we cannot point you to a better solution.

It is better not to use the term "vectorize" as that has a different meaning in the context of R programming; plus you could also "vectorize" raster data to lines or points, so the term is not very clear.

I see. I can't run time(x) <- date. It shows an error saying that 'time<-' can't be applied to an object of class "SpatRaster". But in any case, I think that by subsetting the names of the files I will be able to store them in different monthly/weekly folders by year and transform them to polygons as you did. Thanks! — Berta_94, Oct 17 '22 at 16:09
You probably have another package loaded that hides `terra::time<-`. In that case, you can do `terra::time(x) <- date`. It is important to load as few packages as possible, and to inspect the messages emitted when running `library(apackage)`. — Robert Hijmans, Oct 17 '22 at 16:47

Group raster files by week/month

1 Answers1