0

I have MODIS data which contains 8 days product. Some of month has 3 images and some of month has four imgaes. So how I can build a R scrip to calculate monthly sum from this data

I am using the MODIS MOD17A2HGF GPP data set which comes as a 500x500 meter grid of the 8 day sum of GPP. The first file ends in _001 which means that it starts on January 1st and goes until the 8th. The next file, _009, then covers the 8th to the 16th of January and on and on. enter image description here

Now I want to sum GPP for months（January to December）.

In this folder, it have 920 files (20 years) which every year has 46 files.

In ArcGIS, I know how to sum with raster calculator.But it will do so much repetitive works if I sum every month.

enter image description here

In R, I tried to do it just like this

library(terra)
a <- rast('MOD17A2HGF_GPP_2001_001.tif')
b <- rast('MOD17A2HGF_GPP_2001_009.tif')
c <- rast('MOD17A2HGF_GPP_2001_017.tif')
d <- rast('MOD17A2HGF_GPP_2001_025.tif')
GPP_January <- a+b+c+d
writeRaster(GPP_January,'F:/test/GPP_January_2001.tif')

It will also do so much repetitive works if I sum every month.I don't know how to select the specified files of month in R loop.

I know how to sum annual GPP

library(tidyverse)
library(terra)

GPP2001 <- list.files(pattern = '2001',full.names = T) %>% 
  rast() %>% sum()

But how can I sum monthly GPP per year in R or ArcGIS?

jackywang
  • 75
  • 6

1 Answers1

0

Get all files/layers

library(terra)
ff <- list.files(pattern="\\.tif$", full.names = T) 
r <- rast(ff)

For each file, extract the year and the day number

# ff <- c('MOD17A2HGF_GPP_2001_001.tif', 'MOD17A2HGF_GPP_2001_009.tif')
s <- gsub('MOD17A2HGF_GPP_', "", basename(ff))
s <- gsub('.tif', "", s)
s <- strsplit(s, "_")
s <- do.call(rbind, s)
s
#     [,1]   [,2] 
#[1,] "2001" "001"
#[2,] "2001" "009"

Get the month and combine with year

d <- as.Date(as.numeric(s[,2]), origin = paste(s[,1], "-01-01", sep = "")) - 1
m <- format(d, "%m")
ym <- paste0(s[,1], m)
ym
#[1] "200101" "200101"

Now you can calculate the monthly average like this

x <- tapp(x, ym, mean)

(to be more precise, the average for the dates in each month; that is not exactly the monthly average, but probably good enough).

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Excuse me, is the last code 'x=tapp(r,ym,sum)'? – jackywang May 16 '22 at 20:34
  • I can't describe how excited I am at the moment. I was inspired by a similar question you answered in 2018 and I thought it was perfect, but that question asked about averages and I don't know how to modify the last line of code. I didn't expect the original author's answer not long after my question, which is so exciting! – jackywang May 16 '22 at 20:38
  • I use `mean` (not `sum`) because you stated that you want to compute the "monthly average". – Robert Hijmans May 16 '22 at 21:12
  • I am sorry about that, it is a mistake. I have corrected it in my description. Thanks a lot for your answer and profession! – jackywang May 16 '22 at 21:18
  • sum is also a bit suspect because the number of layers per month will vary – Robert Hijmans May 16 '22 at 21:39