How to subtract two Rasters with discrete values in R?

Question

I'm working with these two rasters:

clc_forest18
class      : RasterLayer 
dimensions : 16000, 12000, 1.92e+08  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 4e+06, 5200000, 1200000, 2800000  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source     : r_tmp_2022-02-03_154451_11816_01443.grd 
names      : layer 
values     : 0, 1  (min, max)

clc_forest00
class      : RasterLayer 
dimensions : 16000, 12000, 1.92e+08  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 4e+06, 5200000, 1200000, 2800000  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source     : r_tmp_2022-02-03_152439_11816_04240.grd 
names      : layer 
values     : 0, 1  (min, max)

I want to get the difference between the clc_forest18 and clc_forest00 in order to see where there are values 1 in 2018 that weren't in 2000. Value 1 is where there is forest cover. In other words I want to see where the forest came back.

I used:

forest_cover_diff <- overlay(clc_forest18, clc_forest00, fun=function(r1, r2){return(r1-r2)})
plot(forest_cover_diff, main="Difference between 2018 and 2000 a.k.a. Where forest came back")

and also:

forest_cover_diff <- clc_forest18 - clc_forest00 
plot(forest_cover_diff, main="Difference between 2018 and 2000 a.k.a. Where forest came back", axes=FALSE) 

but I get a yellow image. I guess it's because im working with discrete values. How can I get the difference then? My goal is to get an image where only new forest pixels (2018) is colored, to see where a rewilding process happened.

Answer 1

Let's make two tiny rasters of 1s and 0s in slightly different patterns:

clc_forest18 = raster(matrix(c(1,1,0,0),2,2))
clc_forest00 = raster(matrix(c(1,0,1,0),2,2))

plot(clc_forest18 - clc_forest00)

subtracting them gives me this:

ie yellow (0) where they are the same, green (+1) where forest18=1 and forest00=0, and white (-1) where its the other way around.

If you are seeing an all yellow plot then either:

• Your data are equal everywhere (eg if you plot the difference of the same raster: plot(clc_forest18 - clc_forest18))
• You're not seeing the difference because the high resolution of your data means R can't plot all the cells. Try:

fdiff = clc_forest18 - clc_forest00
table(values(fdiff))
    
 -1  0  1 
  1  2  1 

Answer 2

A couple ways of evaluating change between two classified rasters are a proportional evaluation or entropy.

Let's create some example data representing forest/nonforest [0,1].

library(raster)
library(spatialEco)

clc_forest18 <- raster(ncols=100, nrows=100)
  clc_forest18[] <- sample(0:1, ncell(clc_forest18), replace=TRUE)
    clc_forest18 <- focal(clc_forest18, matrix(1,5,5), fun=median)  
clc_forest00 <- raster(ncols=100, nrows=100)
  clc_forest00[] <- sample(0:1, ncell(clc_forest00), replace=TRUE)
   clc_forest00 <- focal(clc_forest00, matrix(1,5,5), fun=median)  
plot(stack(clc_forest18, clc_forest00))

Now we can calculate classes representing types of change. You can play with p (percentage threshold) to adjust the sensitivity. The resulting raster values are 1:5 representing: (1) High Decrease, (2) Decrease, (3) Unchanged, (4) Increase, and (5) High Increase.

p=10
rdiff = clc_forest00 - clc_forest18
  dmin=cellStats(rdiff, stat='min', na.rm=TRUE)
    dmax=cellStats(rdiff, stat='max', na.rm=TRUE)
  dpct <- min(c(0.0, dmin)) * min( max(c( 0.0, p)), 100) 
ipct <- max(c(0.0, dmax)) * min( max(c( 0.0, p)), 100) 
rclass <- function(x) { 
  ifelse( x < dpct, 1, 
    ifelse( x >= dpct & x < 0, 2,
      ifelse( x == 0, 3,
    ifelse( x > 0 & x <= ipct, 4,
  ifelse( x > ipct, 5, NA ))))) }
      
forest.chg <- calc(rdiff, fun=rclass)

arg <- list(at=c(1,2,3,4,5), labels=c("High Decrease","Decrease","Unchanged", 
            "Increase", "High Increase"))
cols=c("cyan", "blue", "white", "yellow", "red")
plot(forest.chg, col=cols, axis.args=arg, main="Change")

Another approach is to derive the entropy for pre and post and then take the delta. The direction of the difference will indicate loss or gain of a feature. Ideally, entropy works better for multi-class problems than binomial but, depending on intent can be a robust measure of loss/gain for binomial process.

f18 <- raster.entropy(clc_forest18, d=5, categorical = TRUE)
f00 <- raster.entropy(clc_forest00, d=5, categorical = TRUE)
ent.delta <- f00 - f18
plot(stack(clc_forest18, clc_forest00, ent.delta))