extracted_code.R

library(knitr)
opts_chunk$set(fig.width=7, fig.align='center', fig.height=7, dpi=72)

# load require R packages. Make sure the are installed before loading them
  library("raster")
  library("rgdal")
  library("lubridate")
  library("zoo")
  library("bfast")

# setwd("C:/Data/M4L/FM/BFAST_monitor_V4") 

ndmi_rasterStack <- brick("NDMI_2016_2020_MOOC_final_subset.tif") # NDMI raster stack
Benchmark_forest_mask <- raster("Benchmark_mask_2018_subset.tif") # Forest mask
median_base_image <- brick("S2_multispectral_final_subset.tif")   # multispectral image
datex <- read.csv("image_dates.csv")                       # acquisition date of images
image_acquistion_dates <- as.Date(datex$Date, format("%Y-%m-%d"))
names(ndmi_rasterStack)[c(1:301)] <- as.character(image_acquistion_dates)

plotRGB(median_base_image, 3, 2, 1, stretch='lin')


## > **Question 1**: What can you tell about the study area from this image? What land cover types do you think are within the study area?



plot((subset(ndmi_rasterStack,1:16)/10000), zlim = c(-1, 1))

## > **Question 2**: Some NDMI layers are blank (e.g. 04.11, 05.21, and 06.30). Why is this the case? What could be the reason for changing NDMI values apart from forest disturbances?



plot(Benchmark_forest_mask)



masked_ndmi_rasterStack <- mask(ndmi_rasterStack,Benchmark_forest_mask, maskvalue =NA)



# set the year you want to start monitoring forest disturbance
start_monitoring_year <- 2020

# Get the pixel number on the specified coordinates
Stable_forest_pixel <- cellFromXY(ndmi_rasterStack, c(37.14771, -16.86283))

# Pull out the NDMI time series for a single pixel for stable forest
ndmi_values2 <- as.numeric(ndmi_rasterStack[Stable_forest_pixel])/10000

# replace zeros with NA in case they exist
ndmi_values2[ndmi_values2 ==0]<- NA

# convert ndmi values to time series object
ndmix2 <- bfastts(ndmi_values2, image_acquistion_dates, type = c("irregular"))

# Apply Bfastmonitor 
bfm <- bfastmonitor(data = ndmix2,formula = response ~ harmon, start=c(start_monitoring_year), history = c("all"), order = 3, level=c(0.001,0.001))



plot(bfm, xlab ="Time", ylab ="NDMI" , main ="")



# set the year you want to start monitoring forest disturbance
start_monitoring_year <- 2020

# Get the pixel number on the specified coordinates
deforested_pixel <- cellFromXY(ndmi_rasterStack, c(37.13116, -16.85228))

# Pull out the NDMI time series for a single pixel for stable forest
ndmi_values <- as.numeric(ndmi_rasterStack[deforested_pixel])/10000

# replace zeros with NA in case they exist
ndmi_values[ndmi_values ==0]<- NA

# convert ndmi values to time series object
ndmix <- bfastts(ndmi_values, image_acquistion_dates, type = c("irregular"))

# Apply Bfastmonitor 
bfm1 <- bfastmonitor(data = ndmix,formula = response ~ harmon, start=c(start_monitoring_year), history = c("all"), order = 3, level=c(0.001,0.001))


plot(bfm1, xlab ="Time", ylab ="NDMI" , main ="")


## > **Question 3**: How does this pixel compare to the first. What type of forest disturbance could it have been?


xbfastmonitor <- function(x, timestamps = image_acquistion_dates, start_monitoring_year=2020) {
    ndmix <- bfastts(x, image_acquistion_dates, type = c("irregular"))
    bfm <- bfastmonitor(data = ndmix,formula = response ~ harmon, start=c(start_monitoring_year), history = c("all"), order = 3, level=c(0.001,0.001))
    return(c(breakpoint = bfm$breakpoint, magnitude = bfm$magnitude))
}




rad <- calc(ndmi_rasterStack/10000, xbfastmonitor, start_monitoring_year = 2020, timestamps = image_acquistion_dates, progress = "text")



plot(rad)


## 
## **Question 4**: BFASTmonitor is quite a sensitive algorithm and mapped many disturbances. Based on the spatial and temporal distribution of the disturbances, do you think they correspond to deforestation? Motivate your answer.


## writeRaster(rad[1], filename = "Time_bfastmonitor.tif", overwrite =T)
## writeRaster(rad[2], filename = "Magnitude_bfastmonitor.tif", overwrite =T)
##