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