1. Introduction

This notebook illustrates land cover change in Montes de Maria from 2001 to 2017. Land cover categories are represented by two raster datasets extracted from the Collection 6 MODIS Land Cover MCD12Q1 product.

MCD12Q1 supplies global maps of land cover at annual time steps and 500-m spatial resolution for 2001-present. The product contains 13 Science Data Sets (see Table 1), including 5 legacy classification schemes (IGBP, UMD, LAI, BGC, and PFT) and a new three layer legend based on the Land Cover Classification System (LCCS) from the Food and Agriculture Organization.

More details on this product can be found in the MCD12_User_Guide_V6

# install.packages("diffeR")
# Loading the package
library(raster)
## Loading required package: sp
library(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: reshape2

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

Change metrics are explained in this article: Pontius & 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 MCD12Q1 land cover datasets were obtained previously from Google Earth Engine.

ruta <- "../lc/MODIS/"
files <- list.files(path = ruta, pattern = "tif", all.files = TRUE)
files
## [1] "montesLC2000.tif" "montesLC2017.tif"
(file1 <- paste(ruta,files[1],sep=""))
## [1] "../lc/MODIS/montesLC2000.tif"
lc1 <- raster(paste(ruta,files[1],sep=""))
lc2 <- raster(paste(ruta,files[2],sep=""))
mlc <- stack(lc1,lc2)
names(mlc) <-  c("lc00", "lc17")

Let’s check the landcover stack:

mlc 
## class      : RasterStack 
## dimensions : 206, 218, 44908, 2  (nrow, ncol, ncell, nlayers)
## resolution : 500, 500  (x, y)
## extent     : -8298109, -8189109, 1025555, 1128555  (xmin, xmax, ymin, ymax)
## crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : lc00, lc17 
## min values :    0,    0 
## max values :  255,  255

Let’s load the montes shapefile:

(montes <- shapefile("../MontesdeMaria_WGS84.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 15 
## extent      : -75.70418, -74.77214, 9.22535, 10.14548  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 11
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                        MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,    Shape_Area,      layer,                                                                             path 
## min values  :         13,      13212,     CHALÁN,                              1780,   83.71001417,      2017,    BOLÍVAR,  0.5720244142, 0.00689090159, bol_montes, /Users/ivanlizarazo/Documents/ivan/PR/montes/bol_montes.shp|layername=bol_montes 
## max values  :         70,      70823,   ZAMBRANO, Ordenanza 6 de Octubre 30 de 1968, 1035.07793303,      2017,      SUCRE, 2.08052417308, 0.08525488217, suc_montes, /Users/ivanlizarazo/Documents/ivan/PR/montes/suc_montes.shp|layername=suc_montes

Let’s reproject the shapefile:

montes.rep <- spTransform(montes, crs(mlc))

Now, crop the land cover data:

# Crop land cover data by extent of Montes de Maria
mlc.crop <- crop(mlc, extent(montes.rep))

Now, mask the land cover stack:

mlc.mask <- mask(x = mlc.crop, mask = montes.rep)

Now, 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 : 204, 251, 51204, 2  (nrow, ncol, ncell, nlayers)
## resolution : 500, 500  (x, y)
## extent     : 1146321, 1271821, 1515181, 1617181  (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     : memory
## names      : lc00, lc17 
## min values :    2,    2 
## max values :   17,   17
unique(mlc2[[2]])
## [1]  2  4  8  9 10 11 12 13 17

4. Visualize Land Cover Change

library(leaflet)
library(RColorBrewer)
pal <- colorFactor(c('#152106', '#225129', '#369b47', '#30eb5b', '#387242', 
    '6a2325', '#c3aa69', '#b76031', '#d9903d', '#91af40', 
    '111149', '#cdb33b', '#cc0013', '#33280d', '#d7cdcc',
    '#f7e084', 
    '#aec3d4' ), 
    c(1,  2,  3,  4, 5,  6,  7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17), na.color = "transparent")

Let’s remind what are the IGBP classes:

library(knitr)
include_graphics("../lc/IGBP_classes.png")

LC 2000

leaflet() %>% addTiles() %>%
  addRasterImage(mlc2[[1]], colors = pal, opacity = 1.0, method="ngb") %>%
  addLegend(pal = pal, values = c(1,  2,  3,  4, 5,  6,  7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17),
    title = "Montes de María - LandCover 2000")

LC 2017

leaflet() %>% addTiles() %>%
  addRasterImage(mlc2[[2]], colors = pal, opacity = 1.0, method="ngb") %>%
  addLegend(pal = pal, values =  c(1,  2,  3,  4, 5,  6,  7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17),
    title = "Montes de María - LandCover 2017")

5. Metrics of land cover change from 2000 to 2017

comp <- mlc2[[2]] # 2017 Land Cover
ref <-  mlc2[[1]] # 2000 Land Cover
ctmatCompRef <- crosstabm(comp, ref)
ctmatCompRef
##      2   4    8    9   10  11 12  13   17
## 2  291  23  214   14    4  11  0   0    0
## 4   34 121  104   16    0   0  0   0    0
## 8  304 332 6232 2002   95  10  0   0    0
## 9    5  26  726 9413  500   5  0   0    0
## 10   1   5  103  761 2272   8  1   0    0
## 11   1   0    9    2    1 108  0   0    1
## 12   0   0    0    1    4   0  0   0    0
## 13   0   0    0    3    4   0  0 292    0
## 17   0   0    0    0    0   0  0   0 1607
categoryComponentsPlot(ctmatrix = ctmatCompRef, units = "pixels")

The above figure shows that the most significant changes involved the following classes: * Class 8 - Woody Savannas * Class 9 - Savannas * Class 10 - Grasslands * Class 11 - Permanent Wetlands

It shows that there were also changes in: * Classes 2 and 4 - Broadleaf Forests * Class 12 - Cropland

differenceMR(comp, ref, eval='original', percent=TRUE)
##   Resolution Quantity Exchange     Shift  Overall
## 1        500 7.282007 13.23151 0.2532533 20.76677
diffTablej(ctmatCompRef)
##    Category Omission Agreement Comission Quantity Exchange Shift
## 1         2      345       291       266       79      488    44
## 2         4      386       121       154      232      286    22
## 3         8     1156      6232      2743     1587     2296    16
## 4         9     2799      9413      1262     1537     2498    26
## 5        10      608      2272       879      271     1196    20
## 6        11       34       108        14       20       26     2
## 7        12        1         0         5        4        2     0
## 8        13        0       292         7        7        0     0
## 9        17        1      1607         0        1        0     0
## 10  Overall     5330     20336      5330     1869     3396    65

In order to better understand individual land cover changes, we used the Google Earth Engine (GEE) platform.

These were the results:

library(knitr)
include_graphics("../lc/MODIS_Change_2000_2017.png")

6. Identify hotspots of change

Where are located the most important changes?

The following figures can provide clues:

Large hotspots

Change regions greater than one hundred hectares are shown:

library(knitr)
include_graphics("../lc/pngs1206/polygons3.png")

LC 2001

library(knitr)
include_graphics("../lc/pngs1206/IGBP2001.png")

LC 2017

library(knitr)
include_graphics("../lc/pngs1206/IGBP2017.png")

Polygons with Labels

library(knitr)
include_graphics("../lc/pngs1206/polygons1.png")

Polygons with Image Background

library(knitr)
include_graphics("../lc/pngs1206/polygons2.png")

7. A preliminary list of sites

The sites to study could be defined around the following list of locations:

That’s all. It was the story told by MODIS Land Cover MCD12Q1 - IGBP - at 500 m spatial resolution.