0

I have in my project to analyze a time series of NDVI data for significant changes. For this I did a pixelwise check:

# Loading all necessary libraries
library(sp)
library(raster)
library(rgdal)
library(rastervis)
library(rgl)
library(ggplot2)
library(ggmap)
library(spatial.tools)
library(tidyverse)
library(rgeos)
library(plyr)
library(dbplyr)
library(quantmod)

#Set the path
setwd("/home/r-Daten")

#Create a rasterstack
NDVI_avhrr_path <- "/home/r-Daten"
all_NDVI <- list.files(NDVI_path,full.names = TRUE,pattern = ".tif$")

#fill the rasterstack with the rasters
NDVI_s <- stack(all_NDVI)

#load a single raster for cell size adjustment with the DTM
rastertest <- raster("VHP.P1981040.tif")

#load the DTM
dgm <- raster("/home/r-Daten/srtm_23_16.tif")


#Adjusting the cell size of the DTM to the cell size of rastertest
dgm_2 <- spatial_sync_raster(dgm, rastertest, method = "bilinear", size_only = FALSE,
                         raster_size, verbose = FALSE)

#Cutting the extend DTM into the actual study area
dgm_2_crop <- crop(dgm_2, dgm)


#Cutting the Rasterstack into the extent of dgm_2_crop
NDVI_s_crop <- crop(NDVI_s, dgm_2_crop)

#Compilation of annual totals since autocorrelation with monthly values could be present. For testing purposes, only one image per year is currently used
fun <- function(x){
  NDVI_st_crop.ts = ts(x,start=c(1981), end=c(2015),frequency = 1)
  x<- aggregate(NDVI_s_crop.ts)
}

#Then the slope is calculated to obtain the direction and magnitude of the trends multiplied by the number of years to obtain the change in NDVI units
NDVI_s_crop.sum <- calc (NDVI_s_crop, fun)
NDVI_s_crop.sum=NDVI_s_crop.sum

#Now we have to see which trends are significant on both sides. So first we extract the p-value:
time <- 1:nlayers(NDVI_s_crop.sum) 
fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[2] }}
NDVI_s_crop.slope=calc(NDVI_s_crop.sum, fun)
NDVI_s_crop.slope=NDVI_s_crop.slope*34

NASCS = NDVI_s_crop.slope
plot(NASCS)

#Reading out the values of the significance test
NASCS.values=values(NASCS)

#getting a quick view into the data
NASCS.values

#Now I want to read the critical values and have found a function on the    youtube video: https://www.youtube.com/watch?v=8842FnuX4gc , as long as I use the library (quantmod): 
critical_values = c(qt(NASCS(0,025, 29)),qt(NASCS(0,975, 29)))
critical_values 
#The term "29" refers to degrees of freedom from a Youtube video: 
#now I get an error: 
Fehler in NASCS(0, 25, 29) : could not find function "NASCS"[/code]

But how am I supposed to say r of what I want the respective 0.05 values? Or. how can I classify the significance result in such a way that at the end I get out a plot in which the DGM is shown as a level plot and on the result of the significance test as hatching, where it is just about to have changed significantly does not change significantly.

At the very end, I would like to make statements about the height level, every 1000 meters, in terms of percentage of how much area has changed significantly with the help of the reclassification of the Siknifikanzergebnisse and the reclassified terrain model in altitude levels.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.