Plot times series for stack raster based on a polygons in my area

Question

I am using R language and I would like to create a EVI MODIS time series plot, EVI vs time. I can did method to create a raster stack and plot the average EVI of the total area over time. But I would like to plot for a select polygons : 140 polygons with centroid or lat/long coordinate for centroid for polygon) within the area.

i have 1 raster for one month = 8 raster for year :

These are the polygons link for polygon

https://drive.google.com/drive/u/0/folders/1_frvUnKXF1ASD4c65pHs-asUXJUveyh2

This is my script for plot avrage EVI=meanEVI of total area over times ( EVI/times) :

 rm(list=ls())
library(scales)
library(raster)
library(rgdal)
library(ggplot2)
library(rts)
library(sf)

# Create list of NDVI file paths
all_EVI <- list.files("D:/Rteledetection/Pivots/Doy2000",full.names = TRUE, pattern = ".tif$")
# Create a time series raster stack
EVI_stack <- stack(all_EVI)
# apply scale factor
EVI_EVI_stack <- EVI_stack/10000
EVI_EVI_stack
plot(EVI_EVI_stack)
summary(EVI_EVI_stack)

#Calculate Average NDVI
# calculate mean NDVI for each raster
avg_EVI_stack <- cellStats(EVI_EVI_stack,mean)
# convert output array to data.frame
avg_EVI_stack <- as.data.frame(avg_EVI_stack)

# view column name slot
names(avg_EVI_stack)

#renomer le nom avec meanEVI
names(avg_EVI_stack) <- "meanEVI"

# view cleaned column names
names(avg_EVI_stack)

avg_EVI_stack$product <- "MOD13Q1"

# add a "year" column to our data
avg_EVI_stack$year <- "2000"

julianDays <- gsub(pattern = "MOD13Q1.006__250m_16_days_EVI_doy2000|_aid0001","", row.names(avg_EVI_stack))

head(julianDays)

# add julianDay values as a column in the data frame
avg_EVI_stack$julianDay <- julianDays
# set the origin for the julian date (1 Jan 2000)
origin <- as.Date("2000-01-01")
# convert "julianDay" from class character to integer
avg_EVI_stack$julianDay <- as.integer(avg_EVI_stack$julianDay)
avg_EVI_stack$Date<- origin + (avg_EVI_stack$julianDay-1)

#if run!
head(avg_EVI_stack$Date)
# What are the classes of the two columns now? 
class(avg_EVI_stack$Date)

# plot EVI
f <- ggplot(avg_EVI_stack, aes(x= Date, y=meanEVI), na.rm=TRUE) +
  geom_line(color = "blue", size = 1)+ 
  geom_point(size=2,color="black",pch=1)+
  ggtitle("MODIS EVI - 2000 Mato Grosso  Site") +
  xlab("Date") + ylab("Mean EVI_No units (Scale Factor = 0.0001)")+
  theme(text = element_text(size=10))+theme_classic()+
  scale_x_date(labels = date_format("%m-%Y"))+ 
  scale_y_continuous(name="Mean EVI_No units (Scale Factor = 0.0001)", limits=c(0.3, 0.5))+
  theme_bw()


plot(f)

and these results:

I'm trying to calculate the centroid for my polygons like this :

 #calcul centroid 

    centroid = gCentroid(polyk,byid=TRUE)
    plot(polyk)
    points(coordinates(polyk),pch=1)
    points(centroid,pch=2)

• but is not a good idea i need a simple method to plot times series for my polygon within the area. I know I can us extract::raster but I don't know how and where can I write this line of script !?
• i have idea if i do vector to raster to my polygons and i do mask and i multiply the mask with my raster stack and calculate sum for all raster stack after and plot it but i don't know how i can do this script. Can you help with an explanation and example of what I need?

I know my script is not good but I am trying.

Answer 1

You partitioned your code in a way that obscures what object types you are actually using and what led you to the point that you are stuck. In the future please try to clarify your question before asking, ideally providing example data (note; we do not care about a picture of your polygons). Also, please include the package dependencies, and all of your code that generated your objects. In trying replicating your code I found that in the call to ggplot you nested the date_format function which is from the scales package. We should not have to track this type of thing down in order to provide you help.

Your problem is quite easy to address with the raster::extract function. First, lets create some example data consisting of a stack of 12 [-1 to 1] rasters and two corresponding polygons.

library(sp)
library(raster)

r <- raster(ncol = 100, nrow = 100)
  r[] <- runif(ncell(r), -1, 1)
    r <- stack(r)
  for(i in 1:11){
    x=r[[1]]; x[] <- runif(ncell(x), -1, 1)
    r <- addLayer(r, x)  
  }
  names(r) <- paste0("evi_", 1:12)     
polys <- spPolygons(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)), 
                    rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)))
plot(r[[1]])
  plot(polys, add=TRUE)

Now, we can extract the underlying raster values for each polygon, at each time-step. By passing an argument to the "fun" argument we can collapse the raster values to a summary statistic for each polygon. The data is ordered meaning each row of the resulting matrix corresponds to the corresponding polygon (ie., row 1 in the matrix corresponds to polygon 1).

( r.med <- extract(r, polys, median) )          

You can use ggplot to present the results however, in my example I am using the base plot function to iterate each polygon and plot the timeseries. Note that I am pulling the row from the r.med object.

par(mfrow=c(2,1))
  for(i in 1:nrow(r.med)) {
    plot(seq(as.Date("2018/1/1"), by="month", length.out=12), 
         r.med[i,], type="b", xlab="Median EVI", ylab="Month",
         main=paste0("Polygon - ", i) ) 
}