Land Cover Change in Montes de Maria - ESA CCI
Ivan Lizarazo
11.05.2020
1. Introduction
This notebook illustrates land cover change in Montes de Maria from 2000 to 2018. Land cover categories are represented by two raster datasets extracted from the global land cover (LC) products released by the European Space Agency (ESA) within the climate change initiative (CCI).
ESA LandCover CCI supplies global maps of land cover at annual time steps and 300-m spatial resolution for 1992-2018. The maps uses the classification schema shown below.
# install.packages("diffeR")
# Loading the package
library(raster)## Loading required package: splibrary(diffeR)## Loading required package: rgdal## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-2## Loading required package: ggplot2## Loading required package: reshape22. Illustrative Example
The diffeR package will be used for calculation of metrics of difference for comparing two MCD12Q1 datasets representing categorical variables of land cover using the IGBP classification scheme. The datasets represent land cover categories in 2000 and 2017.
To understand change metrics, it is advised to have a look at this article: Robert Gilmore Pontius Jr & Alà Santacruz (2014) Quantity, exchange, and shift components of difference in a square contingency table, International Journal of Remote Sensing, 35:21, 7543-7554, DOI: 10.1080/2150704X.2014.969814
What are quantity, exchange and shift?
Examples 1 (Pontius & Santacruz, 2014)
Examples 2 (Pontius & Santacruz, 2014)
3. Data pre-processing
The complete collection of ESA CCI LC products was downloaded and pre-processed by A.G.-S. Thanks!
ruta <- "/Users/ivanlizarazo/Documents/ivan/pr/montes/lc/ESA_LC/"
files <- list.files(path = ruta, pattern = "tif", all.files = TRUE)
files## [1] "ESA_LC_MontesdeMaria_WGS84_1992.tif"
## [2] "ESA_LC_MontesdeMaria_WGS84_1993.tif"
## [3] "ESA_LC_MontesdeMaria_WGS84_1994.tif"
## [4] "ESA_LC_MontesdeMaria_WGS84_1995.tif"
## [5] "ESA_LC_MontesdeMaria_WGS84_1996.tif"
## [6] "ESA_LC_MontesdeMaria_WGS84_1997.tif"
## [7] "ESA_LC_MontesdeMaria_WGS84_1998.tif"
## [8] "ESA_LC_MontesdeMaria_WGS84_1999.tif"
## [9] "ESA_LC_MontesdeMaria_WGS84_2000.tif"
## [10] "ESA_LC_MontesdeMaria_WGS84_2001.tif"
## [11] "ESA_LC_MontesdeMaria_WGS84_2002.tif"
## [12] "ESA_LC_MontesdeMaria_WGS84_2003.tif"
## [13] "ESA_LC_MontesdeMaria_WGS84_2004.tif"
## [14] "ESA_LC_MontesdeMaria_WGS84_2005.tif"
## [15] "ESA_LC_MontesdeMaria_WGS84_2006.tif"
## [16] "ESA_LC_MontesdeMaria_WGS84_2007.tif"
## [17] "ESA_LC_MontesdeMaria_WGS84_2008.tif"
## [18] "ESA_LC_MontesdeMaria_WGS84_2009.tif"
## [19] "ESA_LC_MontesdeMaria_WGS84_2010.tif"
## [20] "ESA_LC_MontesdeMaria_WGS84_2011.tif"
## [21] "ESA_LC_MontesdeMaria_WGS84_2012.tif"
## [22] "ESA_LC_MontesdeMaria_WGS84_2013.tif"
## [23] "ESA_LC_MontesdeMaria_WGS84_2014.tif"
## [24] "ESA_LC_MontesdeMaria_WGS84_2015.tif"
## [25] "ESA_LC_MontesdeMaria_WGS84_2016_change.tif"
## [26] "ESA_LC_MontesdeMaria_WGS84_2016.tif"
## [27] "ESA_LC_MontesdeMaria_WGS84_2017_change.tif"
## [28] "ESA_LC_MontesdeMaria_WGS84_2017.tif"
## [29] "ESA_LC_MontesdeMaria_WGS84_2018_change.tif"
## [30] "ESA_LC_MontesdeMaria_WGS84_2018.tif"We will select only five years to study land cover changes:
ruta <- "/Users/ivanlizarazo/Documents/ivan/pr/montes/lc/fiveyears/"
(file1 <- paste(ruta,files[1],sep=""))## [1] "/Users/ivanlizarazo/Documents/ivan/pr/montes/lc/fiveyears/ESA_LC_MontesdeMaria_WGS84_1992.tif"files <- list.files(path = ruta, pattern = "tif", all.files = TRUE)
files## [1] "MLC_2000.tif" "MLC_2005.tif" "MLC_2010.tif" "MLC_2015.tif" "MLC_2018.tif"lc1 <- raster(paste(ruta,files[1],sep=""))
lc2 <- raster(paste(ruta,files[2],sep=""))
lc3 <- raster(paste(ruta,files[3],sep=""))
lc4 <- raster(paste(ruta,files[4],sep=""))
lc5 <- raster(paste(ruta,files[5],sep=""))(mlc <- stack(lc1,lc2,lc3,lc4,lc5))## class : RasterStack
## dimensions : 414, 429, 177606, 5 (nrow, ncol, ncell, nlayers)
## resolution : 0.002777778, 0.002777778 (x, y)
## extent : -75.82778, -74.63611, 9.144444, 10.29444 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : MLC_2000, MLC_2005, MLC_2010, MLC_2015, MLC_2018
## min values : 10, 10, 10, 10, 0
## max values : 210, 210, 210, 210, 255Let’s load the montes shapefile:
montes <- shapefile("/Users/ivanlizarazo/Documents/ivan/pr/montes/MontesdeMaria_WGS84.shp")Let’s clip the land cover data:
# Crop land cover data by extent of Montes de Maria
mlc.crop <- crop(mlc, extent(montes))Now, mask the land cover stack:
mlc.mask <- mask(x = mlc.crop, mask = montes)Plot the masked stack:
plot(mlc.mask)
Let’s reproject the land cover dato into the magna sirgas EPSG 3115:
newproj <- "+proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# we need the rgdal package for this
if (require(rgdal)) {
#simplest approach
#pr1 <- projectRaster(mlc, crs=newproj)
# alternatively also set the resolution
mlc2 <- projectRaster(mlc.mask, crs=newproj, res=500, method='ngb')
}mlc2## class : RasterBrick
## dimensions : 215, 216, 46440, 5 (nrow, ncol, ncell, nlayers)
## resolution : 500, 500 (x, y)
## extent : 1147867, 1255867, 1509495, 1616995 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : /private/var/folders/l_/1nv5l0dx78g_d3k3_t85lrjr0000gn/T/RtmpEqo6w2/raster/r_tmp_2020-05-11_191115_89799_06924.grd
## names : MLC_2000, MLC_2005, MLC_2010, MLC_2015, MLC_2018
## min values : 10, 10, 10, 10, 10
## max values : 210, 210, 210, 210, 210unique(mlc2[[2]])## [1] 10 11 30 40 50 60 100 110 120 130 170 180 190 2104. Visualize Land Cover Change
library(leaflet)
library(RColorBrewer)
pal <- colorFactor(c("#AAF0F0", "#FFFF64", "#DCF064",
"#CDCD66", "#006400", "#00A000", "#AAC800", "#003C00", "#788200", "#8CA000", "#BE9600", "#966400", "#FFB432", "#FFEBAF", "#00785A", "#009678", "#C31400", "#FFF5D7", "#0046C8", "#FFFFFF" ), c(10, 11, 20, 30, 40, 50, 60, 70, 100, 110, 120, 130, 140, 150, 160, 170, 190,200, 210, 220), na.color = "transparent")LC 2000
leaflet() %>% addTiles() %>%
addRasterImage(mlc2[[1]], colors = pal, opacity = 1.0, method="ngb") %>%
addLegend(pal = pal, values = c(10, 11, 20, 30, 40, 50, 60, 70, 100, 110, 120, 130, 140, 150, 160, 170, 190,200, 210, 220),
title = "Montes de MarÃa - LandCover 2000")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NALC 2005
leaflet() %>% addTiles() %>%
addRasterImage(mlc2[[2]], colors = pal, opacity = 1.0, method="ngb") %>%
addLegend(pal = pal, values = c(10, 11, 20, 30, 40, 50, 60, 70, 100, 110, 120, 130, 140, 150, 160, 170, 190,200, 210, 220),
title = "Montes de MarÃa - LandCover 2005")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NALC 2010
leaflet() %>% addTiles() %>%
addRasterImage(mlc2[[3]], colors = pal, opacity = 1.0, method="ngb") %>%
addLegend(pal = pal, values = c(10, 11, 20, 30, 40, 50, 60, 70, 100, 110, 120, 130, 140, 150, 160, 170, 190,200, 210, 220),
title = "Montes de MarÃa - LandCover 2010")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NALC 2015
leaflet() %>% addTiles() %>%
addRasterImage(mlc2[[4]], colors = pal, opacity = 1.0, method="ngb") %>%
addLegend(pal = pal, values = c(10, 11, 20, 30, 40, 50, 60, 70, 100, 110, 120, 130, 140, 150, 160, 170, 190,200, 210, 220),
title = "Montes de MarÃa - LandCover 2015")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NALC 2018
leaflet() %>% addTiles() %>%
addRasterImage(mlc2[[5]], colors = pal, opacity = 1.0, method="ngb") %>%
addLegend(pal = pal, values = c(10, 11, 20, 30, 40, 50, 60, 70, 100, 110, 120, 130, 140, 150, 160, 170, 190,200, 210, 220),
title = "Montes de MarÃa - LandCover 2018")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA5. Calculate metrics of land cover change from 2000 to 2018
comp <- mlc2[[5]] # 2018 Land Cover
ref <- mlc2[[1]] # 2000 Land Cover
ctmatCompRef <- crosstabm(comp, ref)ctmatCompRef## 10 11 30 40 50 60 100 110 120 130 170 180 190 210
## 10 2976 0 0 0 1 0 1 0 0 0 1 0 0 2
## 11 0 1460 0 0 0 8 0 0 0 0 0 0 0 0
## 30 0 0 2737 0 1 0 1 0 0 0 0 0 0 0
## 40 0 0 0 913 2 0 1 0 0 0 2 0 0 0
## 50 4 0 8 2 146 0 0 0 1 0 0 0 0 0
## 60 93 2 491 20 0 15037 0 0 2 0 0 1 0 0
## 100 0 0 2 3 0 1 42 0 0 1 1 0 0 0
## 110 0 0 0 0 0 0 0 21 0 0 0 0 0 0
## 120 0 0 0 0 2 1 0 0 46 0 0 0 0 0
## 130 2 0 5 4 0 0 0 0 0 223 0 0 0 0
## 170 2 0 5 1 0 0 0 0 0 0 568 0 0 11
## 180 0 0 0 0 9 0 0 0 0 0 0 252 0 0
## 190 5 2 17 6 0 1 0 1 0 7 5 10 139 8
## 210 3 0 11 3 32 0 0 0 0 1 2 19 0 686categoryComponentsPlot(ctmatrix = ctmatCompRef, units = "pixels")
differenceMR(comp, ref, eval='original', percent=TRUE)## Resolution Quantity Exchange Shift Overall
## 1 500 2.795996 0.1227323 0.2531354 3.171864diffTablej(ctmatCompRef)## Category Omission Agreement Comission Quantity Exchange Shift
## 1 10 109 2976 5 104 8 2
## 2 11 4 1460 8 4 4 4
## 3 30 539 2737 2 537 4 0
## 4 40 39 913 5 34 8 2
## 5 50 47 146 15 32 10 20
## 6 60 11 15037 609 598 6 16
## 7 100 3 42 8 5 4 2
## 8 110 1 21 0 1 0 0
## 9 120 3 46 3 0 4 2
## 10 130 9 223 11 2 0 18
## 11 170 11 568 19 8 8 14
## 12 180 30 252 9 21 0 18
## 13 190 0 139 62 62 0 0
## 14 210 21 686 71 50 8 34
## 15 Overall 827 25246 827 729 32 66