1

I have monthly .tiff images from 2000-2020 in following way:

200001
200002
200003
.
.
.
202012
(format is yyyymm)

Using the code given below, I am able to calculate monthly averages but I am not able to create line time-series plot for a raster stack. I tried using ggplot but it gives an error that dataframe is to be used and not rasterbrick.

library(raster)

setwd("C:/Users/shc-user/Desktop/Data")
rast <- list.files(path="C:/Users/shcser/Desktop/Data",pattern='.*tif',full.names=TRUE)
s <- stack(rast)
names(s)[1:240]
grp <- rep(1:240, rep(3,240))
grp = substr(names(s),1,6)
grp [1:240]
alb_mean <- stackApply(s, grp, mean)

I need help in extending the code forward to plot line time-series of monthly averages where x-axis will show months and y-axis will show data values. I am writing the script in R. An explaination of the coding steps would be helpful.

1
  • One small note, the use of *tif in list.files is not providing the wildcard that you are thinking. R uses regular expressions so, the correct syntax would be pattern="tif$" indicating a "tif" wildcard occurring at the end of the name. Apr 12, 2023 at 17:44

1 Answer 1

1

Well, the information varies by pixel so, what is your plan to represent the temporal process of the entire raster? You certainly do not expect to create a line graph representing the timeseries of all the pixels?

Commonly, timeseires graphs represent the aggregation of a discrete area (eg., a polygon of interest) and not the entire raster. Collapsing the variation of all the pixels into a single measure, of say the mean, would not be very informative. That said, here is a terra based workflow that will plot a monthly timeseries based on the mean of the entire raster (all pixels), infer at your own risk. The process will be slightly different if extracting values for an area so, if this is the case please clarify your question. I am using terra because it is the replacement, by the same developer, for the raster library which is being depreciated.

First, we create a date vector that represents your timeseires (tried to use your date format) then create some dummy raster data with 1-100 bounds.

library(terra)

# Coerce existing vector to monthly
( d <- paste0("2000",c(paste0("0",1:9),10:12)) )
( months <- as.Date(paste0(d,"01"), "%Y%m%d") )

# or, create monthly sequence using start, end
( months <- seq(as.Date(paste0("200001", "01"), "%Y%m%d"),  
              as.Date(paste0("200312", "01"), "%Y%m%d"), 
              "month") )

suppressWarnings({
  s <- rast(nrows=100, ncols=100, nlyrs=length(months))
    s <- rast(lapply(1:nlyr(s), function(i) {
           mm <- c(sample(1:50,1), sample(51:100,1))   
           s[[i]] <- setValues(s[[i]], runif(size(s[[i]]), 
                               mm[1], mm[2] )) 
       }))
  })
       plot(s[[1:6]]) # plots first 6 rasters in timeseires

Now, we pull the mean for each raster and plot the results. Note, the bracket index is simply pulling the first column in the data.frame resulting from the terra::global function.

( means <- global(s, fun="mean")[,1] )
  plot(months, means, type="b", pch=20)  

Or, using ggplot

library(ggplot2)
dat <- data.frame(date=months, ndvi=means)
ggplot(dat) +
  aes(x = date, y = ndvi) +
  geom_line( ) +
  scale_x_date(date_breaks="4 month", date_labels="%B-%Y")

plot of monthly means

Note; to create the months vector you can use as.Date to coerce your existing data. However, the as.Date needs day in the format so, with only monthly strings you need to perform some small string manipulation (adding "01" to each month) to make the function work.

d <- c("200001", "200002", "200003") 
( d <- paste0(d, "01") )
( months <- as.Date(d,"%Y%m%d") )
7
  • HI @Evans, thanks for the response. I am new to R. Could you explain me the relevance of the following lines. (i) mm <- c(sample(1:50,1), sample(51:100,1)). (ii) What does mm[1], mm[2] mean? (iii) In plot(s[[1:10]]), are we trying to plot 1 to 10 layers? And if I want to plot time series from 2000-2020 and I have 240 files, then, should I write plot(s[[1:240]])? (iv) Could you also please explain the relevance of [,1] and type="b", pch=20 in following lines, respectively( means <- global(s, fun="mean")[,1] ) plot(months, means, type="b", pch=20) Apr 5, 2023 at 15:23
  • 1
    Sorry, a few lines of code for global and plot got un-annotated. Ignore the block of code with sample and the mm object, it is for creating example data and not directly relevant to your question. Take a look at ?par and ?plot to understand relevant arguments eg., type controls plot being point, line or both and pch is the type of point symbol. Using plot(s[[1:10]]) is simply plotting the first 10 rasters (you really do not want to attempt to plot all of them) and is independent of the timeseries plot. Apr 5, 2023 at 18:57
  • @Evans I am able to understand the significance of type, pch. But there are few more things that I am still not able to understand. Also, I want to know something about the results which I am finding strange. It would be helpful if you can explain me a bit more regarding the below mentioned points or modify the code accordingly. Apr 5, 2023 at 23:29
  • (i) I modified the code in this way: s <- rast(lapply(1:nlyr(s), function(i) { s[[i]] <- setValues(s[[i]], runif(size(s[[i]]))) })) After this, it gives me the warning that first raster was empty and was ignored, however, my raster is not empty and it is giving the mean value of first raster also. And, if I run my code from second raster, it still gives me the same warning. Why is it so? (ii) My raster is very small with size 12rows X 16columns, so I want to plot the mean of entire raster and not just the first column. How can I do that? Apr 5, 2023 at 23:29
  • Regarding the results: After modifying the lapply code block, first, I plot the mean values for Jan 2006 to Oct 2006. It gave me certain mean values for each month. Then, I plot the mean values for Feb 2006 to Oct 2006. It again gave me certain mean values for each month but it was different from the former. But how's that possible? The mean value for Feb 2006 should always remain the same. Could you please help me with this. Is there anything that I am not doing correctly. Apr 5, 2023 at 23:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.