R skript k analýze trendů v časových řadách NDVI

# install pacman library if needed
if(!require("pacman")){install.packages("pacman")}
## Loading required package: pacman
# Load the necessary libraries
pacman::p_load(char = c('reticulate', 'bfast', 'rgee', 'remotes', 'sf', 'magrittr', 'tigris', 'tibble', 'stars', 'raster',
                         'dplyr', 'forecast', 'stars', 'st', 'lubridate', 'imputeTS', 'leaflet', 'classInt',
                        'RColorBrewer', 'ggplot2', 'googledrive', 'geojsonio', 'ggpubr', 'readr', 'zoo',
                        'tseries', 'trend', 'Kendall', 'rgdal', 'rgeos', 'sp', 'mapview', 'caret', 'devtools'), 
                        install = T, update = F, character.only = T)
# Read in the CSV file (Use read_csv from readr package
data <- read_csv("C:\\Users\\prochazka5\\Desktop\\olomouc_single_polygon_na_values_deleted.csv")
## Rows: 420 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): SEGMENT_ID, NDVI, t
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Convert the data to a time series object (), rok a den startu měření, frekvence 12 (měsíční, nebo jiná) 
#(je třeba správně nastavit začátek a konec dle datasetu!)
data_ts <- ts(data, start = c(1984, 01), end = c(2022, 12), frequency = 12)

Including Plots

Prostý graf časové řady NDVI vybraného polygonu v Olomouckém kraji

#bfast 01 detekuje pouze jediny zlom v casove rade, je zameren na detekci jedineho vyznamneho zlomu
bfast_single <- bfast01(data_ts [,2])

Detekce jediného zlomu v časové řadě

Zlom časové řady NDVI vybraného polygonu v Olomouckém kraji

plot(bfast_single)

#bfast originalni algoritmus
bfast_original <- bfast(
  data_ts [,2],
  h = 0.15,
  season = c("harmonic"),
  max.iter = 10,
  breaks = NULL,
  hpc = "none",
  level = 0.05,
  decomp = c("stl"),
  type = "OLS-MOSUM"
)

BFAST algoritmus a detekce zlomů v trendech

Originální bfast algoritmus hledá libovolné množství zlomů, několikrát iteruje celý proces hledání.

plot(bfast_original)

## NULL
#extrakce informaci z outputu bfast
output_zprava <- print(bfast_original)
## 
##   TREND BREAKPOINTS
##   Confidence intervals for breakpoints
##   of optimal 2-segment partition: 
## 
## Call:
## confint.breakpointsfull(object = bp.Vt, het.err = FALSE)
## 
## Breakpoints at observation number:
##   2.5 % breakpoints 97.5 %
## 1   240         276    279
## 
## Corresponding to breakdates:
##   2.5 %    breakpoints 97.5 % 
## 1 2003(12) 2006(12)    2007(3)
## 
##   SEASONAL BREAKPOINTS:  None
mesic_zlomu <- bfast_original$Time
print(mesic_zlomu)
## [1] 276

Breakpoint tedy s 97 procentní pravděpodobností leží mezi lety 2003 až 2007

Aktuální stav

Vylepšený Earth Engine skript, který exportuje do CSV hodnoty NDVI za všechny tři landsat družice (dříve byly separátní skripty)

Skript obsahoval neobvykle velké množství NA hodnot (bylo jich několikanásobně více, než jich mělo být - vyřešeno vymazáním)

Vegetační období - skript byl upraven tak, aby obsahoval data pouze za období Duben - Září

Další práce

1 - aplikace skriptu na celý CSV soubor s NDVI daty za roky 1984-2022
2 - “zpravidelnění” časové řady = měsíční mediánové hodnoty
3 - detekce outlierů
4 - validace výsledků??
5 - rekonstrukce kompletní časové řady za pomoci fúze dat MODIS v GEE za posledních několik let (viz články z let 2021 a 2022)