Examine differences between pairs of raster maps for given years (e.g. between classifications)

library(tidyverse)
library(caret)   #for confusionMatrix
library(diffeR)  #for map comparison
library(knitr)
library(raster)
library(rasterVis)
library(gridExtra)


####FUNCTIONS
nat2pas <- function(x,y) { x == 1 & y == 5 }
pas2nat <- function(x,y) { x == 5 & y == 1 }

pas2oag <- function(x,y) { x == 5 & y == 2 }
oag2pas <- function(x,y) { x == 2 & y == 5 }

nat2oag <- function(x,y) { x == 1 & y == 2 }
oag2nat <- function(x,y) { x == 2 & y == 1 }

ag2oag <- function(x,y) { x == 3 & y == 2 }
oag2ag <- function(x,y) { x == 2 & y == 3 }



binRatify <- function(ras){
  
  ras <- ratify(ras)
  rat <- levels(ras)[[1]]
  rat$code <- c("No Change","Change")
  #if(length(rat) == 1) rat$code <- "unknown"
  #if(length(rat) == 2) rat$code <- c("No Change","Change")
  
  levels(ras) <- rat

  return(ras)
  
}



makeObsLUmap <- function(LU, year) {
  
  
  #add categories for later plotting etc. (see https://stackoverflow.com/a/37214431)
  LU <- ratify(LU)     #tell R that the map raster is categorical 
  rat <- levels(LU)[[1]]    #apply the levels (i.e. categories) 

  uLU <- unique(LU) 

  LUcols <- c()
  LUlabs <- c()
  
  if(1 %in% uLU) { 
    LUcols <- c(LUcols, 'forestgreen') 
    LUlabs <- c(LUlabs, 'Nat')  }
  if(2 %in% uLU) { 
    LUcols <- c(LUcols, 'darkcyan') 
    LUlabs <- c(LUlabs, 'OAg') }
  if(3 %in% uLU) { 
    LUcols <- c(LUcols, 'wheat2') 
    LUlabs <- c(LUlabs, 'Ag') }
  if(4 %in% uLU) { 
    LUcols <- c(LUcols, 'black') 
    LUlabs <- c(LUlabs, 'Oth') }
  if(5 %in% uLU) { 
    LUcols <- c(LUcols, 'orange2') 
    LUlabs <- c(LUlabs, 'Pas') }
  
  rat$LandUse <- LUlabs  
  levels(LU) <- rat 
  
  p <- levelplot(LU, att = "LandUse", col.regions=LUcols, main = paste0(year),
      par.settings = list(
        layout.heights = list(top.padding = 0, bottom.padding = 0),
        layout.widths = list(left.padding = 0, right.padding = 0) 
        )
      )  

  return(p)
}
pathA <- "C:/Users/k1076631/Google Drive/Shared/Crafty Telecoupling/Data/LandCover/MapBiomas23/ClassificationComparison/ASCII/"

pathB <- "C:/Users/k1076631/Google Drive/Shared/Crafty Telecoupling/Data/CRAFTYInput/Data/ObservedLCmaps/"


Aname <- "Original"
Bname <- "Planted Data"

data_yrs <- seq(2001, 2015, 1)

for(i in seq_along(data_yrs)){

  print(data_yrs[i])
  #create a stack of rasters for first set of maps
  innameA <- paste0(pathA,"brazillc_",data_yrs[i],"_PastureB.asc")
  lcA <- raster(innameA)
  
  if(i == 1) sA <- stack(lcA)  else  sA <- stack(sA, lcA)
  
  #create a stack of rasters for second set of maps
  innameB <- paste0(pathB,"NewAgri_brazillc_",data_yrs[i],"_PastureB.asc")
  lcB <- raster(innameB)
  
  if(i == 1) sB <- stack(lcB)  else sB <- stack(sB, lcB) 
  
}

#mask to study area
munis.r <- raster("C:/Users/k1076631/Google Drive/Shared/Crafty Telecoupling/Data/CRAFTYInput/Data/sim10_BRmunis_latlon_5km_2018-04-27.asc")

sA <- mask(x=sA, mask=munis.r)   
sA <- trim(sA, padding=2)

sB <- mask(x=sB, mask=munis.r)   
sB <- trim(sB, padding=2)
LCnames <- c("Nat", "OtherAgri", "Agri", "Other", "Pasture")  #used to label error matrix in loop below

comparisons <- c(nat2pas, pas2nat, pas2oag, oag2pas, nat2oag, oag2nat, ag2oag, oag2ag) 
comparisons_n <- c("nat2pas", "pas2nat", "pas2oag", "oag2pas", "nat2oag", "oag2nat", "ag2oag", "oag2ag")

comparisons <- c(nat2pas, pas2nat, pas2oag, oag2pas, nat2oag, oag2nat, ag2oag, oag2ag) 
comparisons_n <- c("nat2pas", "pas2nat", "pas2oag", "oag2pas", "nat2oag", "oag2nat", "ag2oag", "oag2ag")


for(i in seq_along(data_yrs)){
  
  #i <- 2 #for testing
  
  lul <- list()  #this will hold the plots for the LU map for this year
  lul[[1]] <- makeObsLUmap(sA[[i]], Aname)
  lul[[2]] <- makeObsLUmap(sB[[i]], Bname)
  

  #output the crosstab  
  cat("  \n","  \n","Crosstab ",data_yrs[i],"  \n") 
  xtab <- crosstabm(sA[[i]], sB[[i]])
  colnames(xtab) <- LCnames
  rownames(xtab) <- LCnames
  cat("  \n")
  print(kable(xtab))
  cat("  \n")
  
  pl <- list()  #this will hold the plots for the all map for this year
  ts <- stack(sA[[i]], sB[[i]]) #stack used when creating difference maps

  for(j in seq_along(comparisons)){
    
    #j <- 2 #for testing
    np <- raster::overlay(x=ts, fun=comparisons[[j]])
    
    if(cellStats(np, sum) > 0) #if there are any differences
    {
      np <- binRatify(np) 
      mycols <- c("lightgray", "red")
    } else {
      mycols <- c("lightgray") 
    }
    
        #create the plot
        p <- levelplot(np,
        contour=F, 
        margin=F,
        colorkey=F,
        scales=list(draw=FALSE),
        col.regions= mycols,
        main = comparisons_n[j],
        par.settings = list(
          layout.heights = list(top.padding = 0, bottom.padding = 0),
          layout.widths = list(left.padding = 0, right.padding = 0) 
          )
        )
    
      pl[[j]] <- p   
    
  }
  
  
  print(marrangeGrob(lul, nrow = 1, ncol = 2, top = paste0("Observed Land Use")))
  print(marrangeGrob(pl, nrow = 2, ncol = 4, top = paste0(data_yrs[i])))

}

Crosstab 2001

Nat OtherAgri Agri Other Pasture
Nat 61618 0 0 0 0
OtherAgri 0 2147 1926 0 21838
Agri 0 2010 6542 0 0
Other 0 0 0 1918 0
Pasture 0 0 0 0 64131

Crosstab 2002

Nat OtherAgri Agri Other Pasture
Nat 60922 0 0 0 0
OtherAgri 0 2365 1940 0 21286
Agri 0 2133 6994 0 0
Other 0 0 0 2037 0
Pasture 0 0 0 0 64443

Crosstab 2003

Nat OtherAgri Agri Other Pasture
Nat 60258 0 0 0 0
OtherAgri 0 2574 2022 0 20443
Agri 0 2237 7361 0 0
Other 0 0 0 2038 0
Pasture 0 0 0 0 65187

Crosstab 2004

Nat OtherAgri Agri Other Pasture
Nat 59473 0 0 0 0
OtherAgri 0 2670 2084 0 19461
Agri 0 2526 7470 0 0
Other 0 0 0 2000 0
Pasture 0 0 0 0 65741

Crosstab 2005

Nat OtherAgri Agri Other Pasture
Nat 58886 0 0 0 0
OtherAgri 0 2470 2010 0 18897
Agri 0 3047 7621 0 0
Other 0 0 0 2021 0
Pasture 0 0 0 0 66473

Crosstab 2006

Nat OtherAgri Agri Other Pasture
Nat 58718 0 0 0 0
OtherAgri 0 2328 1765 0 17428
Agri 0 2673 8521 0 0
Other 0 0 0 2052 0
Pasture 0 0 0 0 67940

Crosstab 2007

Nat OtherAgri Agri Other Pasture
Nat 58990 0 0 0 0
OtherAgri 0 2428 1609 0 18039
Agri 0 2922 8607 0 0
Other 0 0 0 2072 0
Pasture 0 0 0 0 66758

Crosstab 2008

Nat OtherAgri Agri Other Pasture
Nat 58810 0 0 0 0
OtherAgri 0 2864 1693 0 18590
Agri 0 2939 8567 0 0
Other 0 0 0 2075 0
Pasture 0 0 0 0 65887

Crosstab 2009

Nat OtherAgri Agri Other Pasture
Nat 58732 0 0 0 0
OtherAgri 0 3175 1553 0 19748
Agri 0 2924 8807 0 0
Other 0 0 0 2093 0
Pasture 0 0 0 0 64393

Crosstab 2010

Nat OtherAgri Agri Other Pasture
Nat 58852 0 0 0 0
OtherAgri 0 3081 1453 0 20093
Agri 0 3350 8576 0 0
Other 0 0 0 2120 0
Pasture 0 0 0 0 63900

Crosstab 2011

Nat OtherAgri Agri Other Pasture
Nat 58830 0 0 0 0
OtherAgri 0 3537 1595 0 20660
Agri 0 2745 9032 0 0
Other 0 0 0 2173 0
Pasture 0 0 0 0 62853

Crosstab 2012

Nat OtherAgri Agri Other Pasture
Nat 58709 0 0 0 0
OtherAgri 0 3375 1475 0 19875
Agri 0 2957 9546 0 0
Other 0 0 0 2239 0
Pasture 0 0 0 0 63249

Crosstab 2013

Nat OtherAgri Agri Other Pasture
Nat 59271 0 0 0 0
OtherAgri 0 3164 1279 0 15694
Agri 0 4034 10493 0 0
Other 0 0 0 2269 0
Pasture 0 0 0 0 65221

Crosstab 2014

Nat OtherAgri Agri Other Pasture
Nat 59699 0 0 0 0
OtherAgri 0 3165 1484 0 13036
Agri 0 4776 10001 0 0
Other 0 0 0 2293 0
Pasture 0 0 0 0 66971

Crosstab 2015

Nat OtherAgri Agri Other Pasture
Nat 59500 0 0 0 0
OtherAgri 0 2666 1316 0 11912
Agri 0 5513 10998 0 0
Other 0 0 0 2249 0
Pasture 0 0 0 0 67271