library(raster) library(sp) library(rgdal) library(st) library(sf) library(ncdf4) library(gtools) # create a rasterstack r = list.files( pattern='.tif') r <- mixedsort(r) r <- stack(r) # now use this function to create annual sums if you are working with monthly single files # or use the results (single files) from Code_4 (in this case, skip this part) fun <- function(x) { r.ts = ts(x, start=c(2001,1), end=c(2020,12), frequency=12) x <- aggregate(r.ts) } r.sum <- calc(r, fun) r.sum=r.sum/12 plot(r.sum) r.sum <- r # calculate the slope and multiply it by the number of years, here=20 (2001-2020) time <- 1:nlayers(r) fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[2] }} r.slope=calc(r, fun) r.slope=r.slope*20 plot(r.slope) # extract the p-values fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[8] }} p <- calc(r.sum, fun=fun) plot(p, main="p-value") # write the p-values to raster writeRaster(p, filename="pvalue_annual.tiff", format="GTiff", overwrite=TRUE) # extract p-values < 0.05 to get confidence level of 95% and create a mask m = c(0, 0.05, 1, 0.05, 1, 0) rclmat = matrix(m, ncol=3, byrow=TRUE) p.mask = reclassify(p, rclmat) fun=function(x) { x[x<1] <- NA; return(x)} p.mask.NA = calc(p.mask, fun) # mask all insignificant NDVI levels and plot NDVI changes with 95% significance level trend.sig = mask(r.slope, p.mask.NA) plot(trend.sig, main="trend with 95% confidence level") # export to significance raster writeRaster(trend.sig, filename="trend_95perc_confidence.tiff", format="GTiff", overwrite=TRUE) # have a look at the results plot(trend.sig) plot(readOGR(dsn = ".", layer = "europe_coast_4326"), add=TRUE)