I'm working with these two rasters:

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)

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.

  • Have you tried this on some smaller rasters? Make some 10x10 matrices with ones and zeroes and test you get out what you expect. Your yellow image might be because you can't display 16000x12000 pixels and are only seeing a subsample? What values are you actually getting?
    – Spacedman
    Commented Feb 3, 2022 at 18:15

2 Answers 2


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))

enter image description here

plot(clc_forest18 - clc_forest00)

subtracting them gives me this:

enter image description here

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
 -1  0  1 
  1  2  1 
  • Thank you, it helped a lot. I decided to subdivide the image in 3 parts to see the coloured pixels better. I set the colorramppalette as folowing: colorRampPalette(c("red","white","blue"))(100). My problem now is that I have an image almost completely white with some blue and red pixels scattered without any kind of geographical reference point..Should I overlap a new map with just the europe borders? Commented Feb 4, 2022 at 12:36
  • If you coloured it with "grey" in the middle would you see the land shape better? I'm assuming your sea areas are NA which is showing as white as well as the non-changed land...
    – Spacedman
    Commented Feb 4, 2022 at 16:34
  • that's exactly the problem sir. Commented Feb 4, 2022 at 16:52

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].


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.

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))

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.