0

I have two Land Cover (LC) classification raster over the same area and the same cell size. One from the year 2018 and the other from the year 2014. Both raster are classified with the values 1 - 6, where each value represents a different LC type.

Now I would like to calculate the change of LC between 2014 and 2018.

For instant I would like to know how many pixels with LC type 1 in 2014 still have LC type 1, and how many have new LC type 2 or LC type 3 and so on...

The same for pixels with LC type 2, 3, 4, 5 and 6 in 2014.

The output would probably be something like a correspondence matrix.

So far I have looked into the Package ‘lulcc’ and Package ‘diffeR’ but couldn't figure it out yet.

2
  • If you want to visualise this differently, "alluvial" diagrams are great for this. cran.r-project.org/web/packages/alluvial/vignettes/… - if you try and then get stuck ask a new question. – Spacedman Feb 19 '19 at 14:39
  • Thanks for your suggestion. I represent my results now as a "alluvial" diagrams. Very nice! – taeffer Feb 21 '19 at 9:45
3

An easy approach is to compute a cross-table. Give this example:

library(raster)

r <- raster()

set.seed(123)

lc1 <- setValues(r, sample(1:6, 64800, replace = T))
lc2 <- setValues(r, sample(1:6, 64800, replace = T))

s <- stack(lc1,lc2)

rasterVis::levelplot(s)

enter image description here

You just use table() function converting RasterStack to data.frame:

table(as.data.frame(s))
##       layer.2
##layer.1    1    2    3    4    5    6
##      1 1826 1790 1788 1814 1820 1872
##      2 1858 1814 1803 1794 1745 1813
##      3 1790 1878 1881 1843 1787 1823
##      4 1750 1752 1717 1794 1810 1829
##      5 1782 1740 1856 1745 1786 1827
##      6 1775 1778 1743 1775 1838 1764

The result is a table with the number of cross-observations per feature class. To obtain the area per cell, just multiply with pixel size (if you are using a meter or inch CRS)

lcArea <- as.matrix(table(as.data.frame(s)))*res(s)[1]*res(s)[1] # in this case is the same table, pixel resolution is 1x1
0
0

You can use the function crosstable in the package raster.

Assuming that you have the maps 2014, 2018 saved as raster landscape_2014 and landscape_2018, the code would be like that:

contingencyTable <- raster::crosstab(landscape_2014, landscape_2018, long = TRUE)

The argument long define if the results will be in data.frame or table format.

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