Example 1

Example 1

Vegetation index: NDMI

History period: 2000-2010

option: “all”

stack subset: no

Monitoring approach: 2010-2016

option: Full monitoring period approach

Regression model: harmonic order 1

Trend: no

source("~/wur_bfast_workshop/R-scripts/tutorial_0.R")
source("~/wur_bfast_workshop/R-scripts/accuracy_assessment.R")
example_title <- 1
results_directory <- file.path(results_directory,paste0("example_",example_title))
dir.create(results_directory)
log_filename <- file.path(results_directory, paste0(format(Sys.time(), "%Y-%m-%d-%H-%M-%S"), "_example_", example_title, ".log"))
start_time <- format(Sys.time(), "%Y/%m/%d %H:%M:%S")

result <- file.path(results_directory, paste0("example_", example_title, ".grd"))
time <- system.time(bfmSpatial(ndmiStack, start = c(2010, 1),
                               formula = response ~ harmon,
                               order = 1, history = "all",
                               filename = result,
                               mc.cores = detectCores()))

write(paste0("This process started on ", start_time,
             " and ended on ",format(Sys.time(),"%Y/%m/%d %H:%M:%S"),
             " for a total time of ", time[[3]]/60," minutes"), log_filename, append=TRUE)

Post Processing

bfm_ndmi <- brick(result)
#### Change
change <- raster(bfm_ndmi,1)
plot(change, col=rainbow(7),breaks=c(2010:2016))

#### Magnitude
magnitude <- raster(bfm_ndmi,2)
magn_bkp <- magnitude
magn_bkp[is.na(change)] <- NA
plot(magn_bkp,breaks=c(-5:5*1000),col=rainbow(length(c(-5:5*1000))))

plot(magnitude, breaks=c(-5:5*1000),col=rainbow(length(c(-5:5*1000))))

#### Error
error <- raster(bfm_ndmi,3)
plot(error)

#### Detect deforestation
def_ndmi <- magn_bkp
def_ndmi[def_ndmi>0]=NA
plot(def_ndmi)

plot(def_ndmi,col="black", main="NDMI_deforestation")

writeRaster(def_ndmi,filename = file.path(results_directory,paste0("example_",example_title,"_deforestation_magnitude.grd")),overwrite=TRUE)

## class       : RasterLayer 
## dimensions  : 334, 334, 111556  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 374055, 384075, -676185, -666165  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/andrei/wur_bfast_workshop/data/Bfast_results/example_1/example_1_deforestation_magnitude.grd 
## names       : layer.2 
## values      : -3706.921, -0.079172  (min, max)

def_years <- change
def_years[is.na(def_ndmi)]=NA

years <- c(2010,2011,2012,2013,2014,2015,2016,2017)
plot(def_years, col=rainbow(length(years)),breaks=years, main="Detecting deforestation after 2010")

writeRaster(def_years,filename = file.path(results_directory,paste0("example_",example_title,"_deforestation_dates.grd")),overwrite=TRUE)

## class       : RasterLayer 
## dimensions  : 334, 334, 111556  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 374055, 384075, -676185, -666165  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/andrei/wur_bfast_workshop/data/Bfast_results/example_1/example_1_deforestation_dates.grd 
## names       : layer.1 
## values      : 2010.342, 2016.737  (min, max)

Accuracy Asessment

Forest_mask <- raster(file.path(workshop_folder,"data/Fmask_2010_Peru.tif"))
validation_forest_map <- raster(file.path(workshop_folder,"data/Validation_forest_2016.tif"))

The accuracy assessment is done following the paper of Olofsson et al. (2014) “Good practices for estimating area and assessing accuracy of land change”.

sample_size <- calcSampleSize(def_years,Forest_mask,c(0.9,0.7),0.01)
samples <- extractRandomSamples(def_years,Forest_mask,sample_size,results_directory,"samples")
val_sample <- extractValidationValues(validation_forest_map, samples, Forest_mask)
conf_matrix <- assessAcuracy(samples,val_sample)
conf_matrix

##        0      1          
## 0 744.00  31.00  775 0.96
## 1 122.00 573.00  695 0.82
##   866.00 604.00 1470   NA
##     0.86   0.95   NA 0.90