Example 7

Vegetation index: NDMI

History period: 2005-2010

option: “all”

stack subset: yes (from 2005)

Monitoring approach: 2010-2016

option: Sequential 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 <- 7
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"))
subsetTimeStack <- function(timestack,year){
  timespan <- which(timestack@z$time>as.Date(paste0(year,"-01-01"),"%Y-%m-%d"))
  result <- brick(timestack[[timespan]])
  result@z$time <- timestack@z$time[timespan]
  result
}

ndmiStack_example_7 <- subsetTimeStack(ndmiStack,2005)

bfmSpatialSq <- function(start, end, timeStack, outdir, ...){
  bfm_seq <- lapply(start:end,
                    function(year){
                      outfl <- paste0(outdir, "/bfm_NDMI_", year, ".grd")
                      bfm_year <- bfmSpatial(timeStack, start = c(year, 1), monend = c(year + 1, 1), formula = response~harmon,
                                             order = 1, history = "all", filename = outfl, ...)
                      outfl
                    })
}
time <- system.time(bfmSpatialSq(2010,2016,ndmiStack_example_7,results_directory, mc.cores = detectCores()))

calcDefSeqYears2 <- function(outdir,outfile,start,end,parameter_value){
  bfast_result_fnames <- list.files(outdir, pattern=glob2rx('*.grd'), full.names=TRUE)
  yearly_def <- lapply(bfast_result_fnames,function(file_name){
    bfm_year <- brick(file_name)
    bfm_year[[1]][bfm_year[[2]]>0] <- NA
    bfm_year[[2]][is.na(bfm_year[[1]])] <- NA
    bfm_year
  })
  bfm_summary <- yearly_def[[length(yearly_def)]]
  for (i in (length(yearly_def)-1):1) {
    bfm_summary[!is.na(yearly_def[[i]][[1]])] <- yearly_def[[i]][!is.na(yearly_def[[i]][[1]])]
  }
  writeRaster(bfm_summary,file.path(outfile))
}

def_years_2005 <- calcDefSeqYears2(results_directory,result,2010,2016,2005)

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_7/example_7_deforestation_magnitude.grd 
## names       : layer.2 
## values      : -7242.039, -0.01936926  (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_7/example_7_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 827.00  39.00  866 0.95
## 1  77.00 579.00  656 0.88
##   904.00 618.00 1522   NA
##     0.91   0.94   NA 0.92

Example 7

Sabina Rosca

Example 7

Post Processing

Accuracy Asessment