Thursday 7 February 2019

Creating Land Cover Change Classification in R?



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.



Answer



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

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...