# Creating Land Cover Change Classification in R? [closed]

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.

• 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

## 2 Answers

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)
`````` 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)*res(s) # in this case is the same table, pixel resolution is 1x1
``````

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.