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 }
pas2ag <- function(x,y) { x == 5 & y == 3 }
ag2pas <- function(x,y) { x == 3 & y == 5 }
nat2ag <- function(x,y) { x == 1 & y == 3 }
ag2nat <- function(x,y) { x == 3 & y == 1 }
nat2oag <- function(x,y) { x == 1 & y == 2 }
oag2nat <- function(x,y) { x == 2 & y == 1 }
pas2oag <- function(x,y) { x == 5 & y == 2 }
oag2pas <- function(x,y) { x == 2 & y == 5 }
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)
}
path <- "E:/OneDrive - King's College London/Research/Projects/Belmont/Modelling/CRAFTY/Crafty Telecoupling/Data/LandCover/MapBiomas4/BrazilInputMaps/Data/Classified/"
data_yrs <- seq(2001, 2018, 1)
for(i in 1:(length(data_yrs)-1)){
print(paste0(data_yrs[i],"-",data_yrs[i+1]))
#create a stack of rasters for first set of maps
innameA <- paste0(path,"LandCover",data_yrs[i],"_PastureB_Disagg.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(path,"LandCover",data_yrs[i+1],"_PastureB_Disagg.asc")
lcB <- raster(innameB)
if(i == 1) sB <- stack(lcB) else sB <- stack(sB, lcB)
}
#mask to study area
munis.r <- raster("E:/OneDrive - King's College London/Research/Projects/Belmont/Modelling/CRAFTY/Crafty Telecoupling/Data/CRAFTYInput/Data/sim10_BRmunis_latlon_5km.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, nat2ag, ag2nat, pas2oag, oag2pas, nat2oag, oag2nat)
comparisons_n <- c("nat2pas", "pas2nat", "nat2ag", "ag2nat", "pas2oag", "oag2pas", "nat2oag", "oag2nat")
for(i in 1:(length(data_yrs)-1)){
#i <- 2 #for testing
lul <- list() #this will hold the plots for the LU map for this year
lul[[1]] <- makeObsLUmap(sA[[i]], paste0(data_yrs[i]))
lul[[2]] <- makeObsLUmap(sB[[i]], paste0(data_yrs[i+1]))
#output the crosstab
cat(" \n"," \n","Crosstab ",data_yrs[i],"-",data_yrs[i+1]," \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(data_yrs[i],"-",data_yrs[i+1])))
print(marrangeGrob(pl, nrow = 2, ncol = 4, top = paste0("Differences")))
}
Crosstab 2001 - 2002
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 61514 | 172 | 203 | 103 | 2660 |
OtherAgri | 48 | 3657 | 481 | 5 | 110 |
Agri | 49 | 545 | 9679 | 4 | 159 |
Other | 56 | 15 | 13 | 2517 | 81 |
Pasture | 2315 | 277 | 479 | 174 | 77685 |
Crosstab 2002 - 2003
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 60377 | 199 | 216 | 119 | 3071 |
OtherAgri | 57 | 3926 | 598 | 8 | 77 |
Agri | 58 | 393 | 10269 | 8 | 127 |
Other | 63 | 10 | 12 | 2628 | 90 |
Pasture | 2352 | 467 | 838 | 161 | 76877 |
Crosstab 2003 - 2004
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 59604 | 212 | 307 | 70 | 2714 |
OtherAgri | 80 | 4205 | 571 | 5 | 134 |
Agri | 54 | 583 | 11141 | 9 | 146 |
Other | 90 | 26 | 29 | 2658 | 121 |
Pasture | 2518 | 398 | 673 | 142 | 76511 |
Crosstab 2004 - 2005
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 58986 | 192 | 235 | 70 | 2863 |
OtherAgri | 89 | 4323 | 864 | 10 | 138 |
Agri | 92 | 462 | 11945 | 10 | 212 |
Other | 81 | 12 | 13 | 2704 | 74 |
Pasture | 2193 | 332 | 469 | 79 | 76553 |
Crosstab 2005 - 2006
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 58475 | 159 | 192 | 89 | 2526 |
OtherAgri | 123 | 4236 | 795 | 3 | 164 |
Agri | 135 | 642 | 12413 | 8 | 328 |
Other | 64 | 12 | 14 | 2709 | 74 |
Pasture | 2386 | 314 | 366 | 87 | 76687 |
Crosstab 2006 - 2007
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 58377 | 175 | 190 | 92 | 2349 |
OtherAgri | 62 | 4500 | 634 | 7 | 160 |
Agri | 132 | 548 | 12803 | 13 | 284 |
Other | 59 | 11 | 14 | 2771 | 41 |
Pasture | 2405 | 372 | 460 | 113 | 76429 |
Crosstab 2007 - 2008
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 57985 | 204 | 271 | 65 | 2510 |
OtherAgri | 53 | 4720 | 679 | 7 | 147 |
Agri | 125 | 703 | 13007 | 11 | 255 |
Other | 78 | 8 | 19 | 2803 | 88 |
Pasture | 2353 | 397 | 528 | 93 | 75892 |
Crosstab 2008 - 2009
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 57916 | 162 | 210 | 53 | 2253 |
OtherAgri | 93 | 5209 | 515 | 10 | 205 |
Agri | 115 | 688 | 13391 | 15 | 295 |
Other | 93 | 6 | 12 | 2793 | 75 |
Pasture | 2451 | 344 | 338 | 70 | 75689 |
Crosstab 2009 - 2010
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 58043 | 149 | 198 | 67 | 2211 |
OtherAgri | 71 | 5361 | 802 | 6 | 169 |
Agri | 115 | 532 | 13556 | 5 | 258 |
Other | 43 | 12 | 17 | 2806 | 63 |
Pasture | 2212 | 299 | 415 | 102 | 75489 |
Crosstab 2010 - 2011
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 57717 | 154 | 229 | 130 | 2254 |
OtherAgri | 60 | 5629 | 534 | 6 | 124 |
Agri | 108 | 724 | 13957 | 10 | 189 |
Other | 69 | 11 | 13 | 2835 | 58 |
Pasture | 2389 | 389 | 583 | 148 | 74681 |
Crosstab 2011 - 2012
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 57418 | 139 | 237 | 83 | 2466 |
OtherAgri | 93 | 5998 | 636 | 10 | 170 |
Agri | 135 | 941 | 14013 | 7 | 220 |
Other | 109 | 8 | 21 | 2894 | 97 |
Pasture | 2356 | 484 | 691 | 92 | 73683 |
Crosstab 2012 - 2013
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 57107 | 167 | 272 | 84 | 2481 |
OtherAgri | 78 | 6472 | 873 | 6 | 141 |
Agri | 107 | 538 | 14770 | 5 | 178 |
Other | 73 | 15 | 12 | 2873 | 113 |
Pasture | 2292 | 574 | 809 | 96 | 72865 |
Crosstab 2013 - 2014
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 56850 | 183 | 264 | 85 | 2275 |
OtherAgri | 85 | 6598 | 903 | 7 | 173 |
Agri | 127 | 647 | 15669 | 14 | 279 |
Other | 47 | 19 | 17 | 2882 | 99 |
Pasture | 2244 | 554 | 685 | 71 | 72224 |
Crosstab 2014 - 2015
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 56528 | 202 | 180 | 45 | 2398 |
OtherAgri | 94 | 6836 | 918 | 8 | 145 |
Agri | 144 | 778 | 16389 | 6 | 221 |
Other | 97 | 11 | 18 | 2856 | 77 |
Pasture | 2398 | 444 | 569 | 89 | 71550 |
Crosstab 2015 - 2016
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 56625 | 221 | 253 | 53 | 2109 |
OtherAgri | 130 | 7130 | 786 | 13 | 212 |
Agri | 160 | 970 | 16621 | 18 | 305 |
Other | 67 | 6 | 10 | 2853 | 68 |
Pasture | 2697 | 350 | 488 | 87 | 70769 |
Crosstab 2016 - 2017
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 56603 | 223 | 244 | 88 | 2521 |
OtherAgri | 79 | 7574 | 839 | 15 | 170 |
Agri | 117 | 845 | 16930 | 9 | 257 |
Other | 52 | 12 | 19 | 2851 | 90 |
Pasture | 2229 | 307 | 414 | 76 | 70437 |
Crosstab 2017 - 2018
Nat | OtherAgri | Agri | Other | Pasture | |
---|---|---|---|---|---|
Nat | 56173 | 154 | 228 | 98 | 2427 |
OtherAgri | 94 | 7798 | 889 | 10 | 170 |
Agri | 159 | 675 | 17324 | 11 | 277 |
Other | 89 | 19 | 23 | 2822 | 86 |
Pasture | 2478 | 340 | 476 | 99 | 70082 |