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 :

enter image description here

These are the polygons link for polygon


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


# 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

#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

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

# view cleaned column names

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))


# 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!
# What are the classes of the two columns now? 

# plot EVI
f <- ggplot(avg_EVI_stack, aes(x= Date, y=meanEVI), na.rm=TRUE) +
  geom_line(color = "blue", size = 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))+


and these results: enter image description here

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

 #calcul centroid 

    centroid = gCentroid(polyk,byid=TRUE)
  • 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.

  • 1
    What is avg_EVI_stack? Does extract(avg_EVI_stack, centroid) do anything?
    – Spacedman
    Dec 13, 2018 at 15:34
  • averge EVI stack before stack EVI i'm callcul the average for all stack EVI=avg_EVI_stack !! Dec 13, 2018 at 16:41
  • 1
    There are lots of tutorials on extracting raster stack values from polygons, and there have been questions here recently on exactly that, looking a lot like this question. Please try again, and then if you get stuck ask a specific question giving more information about your data (eg what does summary print for it) and error messages.
    – Spacedman
    Dec 13, 2018 at 16:58

1 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.


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(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.

  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) ) 
  • i will edit my question with complet code my problem it i have for one year a polygons so for 6 years a have 6 shapfile and in one shapfil i have 100 polygons so my raster is the modis evi TIF 6 raster for one years , i will edit my question for complet script and information and package Dec 13, 2018 at 22:10
  • im edited my question Dec 13, 2018 at 23:36

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.