The R-package ecochange can assist users in monitoring ecosystem changes by computing biodiversity indicators related to structural essential biodiversity variables, including ecosystem areas, fractal dimension, and information-theory indices. The package processes polygons —defined here as Regions of Interest (ROI)— and Ecosystem Remote Sensing Products (ERSP). ecochange harmonizes ERSP using the geometry of the ROI. The users can provide predefined polygons or implement the in-package functions to download ROI representing Geographic Administrative Data Maps (GADM). The package can also download up to 16 ERSP available in three global databases: Forest Cover (GFC; Hansen et al., 2013), Tree-canopy Cover (TC; Sexton et al., 2013), and Water Surface (GWS; Pekel et al., 2017). Routines of the package integrate these and other predefined rastersif these data cover the extent of the processed ROI. In-package functions for statistical analysis and graphics provide tools to inform and visualize ecosystem changes.


The first requirement is to install and load the package in the R environment. Users can install the package and the dependencies running install.packages(). They can load it using either require() or library(). The package relies on several dependencies. It uses these to improve execution performance during spatial analysis. The users must be sure that all these are correctly installed.

# install.packages('ecochange')

dependencies <- c('ecochange','raster','rgdal','parallel', 'R.utils', 'rvest','xml2',
            'gdalUtilities','rgeos','viridis', 'rasterVis','rlang', 'rasterDT')

# Values in the next list must be allx TRUE
sapply(dependencies, require, character.only = TRUE)

Retrieving spatial data

Function getrsp can download the ERSP by processing the extent of a predefined region of interest (polygon geometry or GADM unit).

Here we implement getrsp() to download a polygon corresponding with a Colombian municipality. We control the function using three arguments: level, country, and mc.cores. These specify the GADM level (level 2 retrieves municipalities), the country ISO, and the number of cores, respectively.

Harmonizing spatial data

Function rsp2ebv processes ERSP to produce harmonized raster-data sections with the cell values enclosed in an ROI.


ebv <- rsp2ebv('Chimichagua', lyrs = c('treecover2000','lossyear'), mc.cores = 2)
## [1] "The new data will be stored in /tmp/RtmpZ8Zxu5/ecochange:"
## class      : RasterStack 
## dimensions : 2149, 2098, 4508602, 2  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 609992.8, 672932.8, 987047.8, 1051518  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs 
## names      : treecover2000, lossyear 
## min values :             0,        0 
## max values :           255,      255

plot of chunk unnamed-chunk-2

Deriving ecological changes

Function echanges can generate ecosystem-change maps by masking cell values in a layer of ecosystem changes over a target set of ecosystem variables.

This function can process the raster-data section previously harmonized with rsp2ebv(). It can also inherit arguments in rsp2ebv() and derive the data sections as an intermediate step.

By default, echanges() process the raster-data sections considering that the first layer is the target variable and the last layer is the change raster. If these layer positions are incorrect, the users must provide two additional arguments: eco and change. Such arguments provide regular expressions matching the names of the variables in the raster-data sections. For instance, we can use the arguments eco and change to provide the expressions 'tree' and 'loss' corresponding with the layers 'treecover2000' and 'lossyear', respectively.

Arguments eco_range and change_vals control the pixel values processed in the target variable and the change raster. Here we maintain the default in eco_range to process canopy covers between 1 and 100% and modify the argument change_vals to process pixels 0, 10, and 19 in the layer 'lossyear'. These last values represent deforested areas observed during the years 2000, 2010, and 2019.

The output of echanges() is a RasterStack object with dimensionality defined by the values provided in change_vals. Here we obtain three layers corresponding with deforestations observed during the three years.

ech <- echanges(ebv, eco = 'tree', echanges = 'loss',
                change_vals = c(0,10,19), mc.cores = 2)

plot of chunk unnamed-chunk-3 The plot suggests that canopy covers between 1-100% has decreased between 2000-2019. Let's compute the areas occupied by each canopy-cover value.

Calculating biodiversity indicators

Gauging indicators

In-package gaugeIndicator() calculates the biodiversity indicators by processing ecosystem-change maps from echanges(). It can also inherit arguments in rsp2ebv() and echanges() to derive the ecosystem-change maps.

Argument metric in the function provides the name of a biodiversity indicator. By default, metric = 'area_ha' makes the function to compute ecosystem areas (hectare). Besides, the users can provide hundreds of other abbreviations available in the landscape dependency (Hesselbarth et al., 2019). Examples of these include 'condent' and 'lsm_l_frac_mn' which make the function to compute conditional entropy (dimensionless) and mean fractal dimension index for all patches in the raster (dimensionless).

The function output is a tibble of indicators ordered according to classes in the ecological variable. This output can be visualized using the in-package plotind().


# computating ecosystem areas (default)
gi <- gaugeIndicator(ech, mc.cores = 2)
## # A tibble: 300 x 6
##    layer level class id    metric  value
##    <chr> <chr> <dbl> <lgl> <chr>   <dbl>
##  1 eco_0 class     1 NA    area_ha 1071.
##  2 eco_0 class     2 NA    area_ha  726.
##  3 eco_0 class     3 NA    area_ha  368.
##  4 eco_0 class     4 NA    area_ha  576.
##  5 eco_0 class     5 NA    area_ha  294.
##  6 eco_0 class     6 NA    area_ha  240.
##  7 eco_0 class     7 NA    area_ha  192.
##  8 eco_0 class     8 NA    area_ha  142.
##  9 eco_0 class     9 NA    area_ha  141.
## 10 eco_0 class    10 NA    area_ha  143.
## # … with 290 more rows

plot of chunk unnamed-chunk-4

Forest areas in the municipality have decreased from to between 2000-2019. We can test if the observed deforestation has affected the landscape diversity by computing changes in conditional entropy.

Sampling indicators across grids

In-package sampleIndicator() divides into fixed-size grids each of the scenes of a stack of ecosystem-spatial data and samples a biodiversity indicator by every grid. This function helps the users to visualize changes in the spatial distribution of the biodiversity indicators.

Default metric = 'condent' makes the function to compute conditional entropy (dimensionless). Users can define any other biodiversity indicator. Here we calculate changes in conditional entropy to measure ecosystem degradation in terms of changes of pixel adjacency and canopy-cover diversity. Default side (missing) can find an optimum grid size (m) that returns at least a numeric value in one grid by scene. The users can change the side default to specify a customized grid size. The function output is a RasterStack object with the same number of scenes of the ecosystem-change map. For instance,

si_ent <- sampleIndicator(ech, mc.cores = 2)

plot of chunk unnamed-chunk-5

Changes in conditional entropy across the ecosystem-change map between 2000-2019 remained relatively constant. Let's compute means and standard deviations.

Calculating statistics

Function EBVstats can map statistics over scenes in RasterStack objects. Here we implement EBVstats() to derive statistics for the changes of entropy observed across the grid maps.

##Statistics for changes in entropy:
sts <- EBVstats(si_ent)
## # A tibble: 3 x 7
##   layer  n.grids   min  mean   max    sd   skew
##   <chr>    <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 eco_0      576  0.44  1.34  1.83 0.349 -0.781
## 2 eco_10     576  0.44  1.35  1.84 0.356 -0.793
## 3 eco_19     576  0.45  1.38  1.89 0.361 -0.797
## In-package bar plot method:

plot of chunk unnamed-chunk-6

Averages and standard deviations for the grids suggest the relationship between changes in pixel adjacency and canopy-cover diversity between 2000-2019 has remained relatively constant across the municipality.


This is a brief introduction to the package ecochange. There are other arguments to control the package functionality. For instance, the function listGP can print detailed information about ERSP downloadable with the package. Likewise, the function echanges() includes other arguments, e.g., 'sp_dist', 'get_unaffected', and 'binary_output'; these arguments help the users to provide an alternative raster indicating a species distribution range, or to switch analyses focussing on affected/unaffected areas, or to format the outputs into binary layers, respectively. Please refer to the package documentation for more information.

Citing ecochange

The package can be cited using the citation() function. Please cite ecochange appropriately in your work.



Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova, S. A. A., Tyukavina, A., … & Kommareddy, A. (2013). High-resolution global maps of 21st-century forest cover change. science, 342(6160), 850-853.

Hesselbarth, M.H., Sciaini, M., With, K.A., Wiegand, K., & Nowosad, J., (2019). landscapemetrics: an open-source R tool to calculate landscape metrics. Ecography.

Pekel, J. F., Belward, A., & Gorelick, N. (2017). Global Water Surface Dynamics: Toward a Near Real Time Monitoring Using Landsat and Sentinel Data. In AGU Fall Meeting Abstracts.

Sexton, J. O., Song, X. P., Feng, M., Noojipady, P., Anand, A., Huang, C., … & Townshend, J. R. (2013). Global, 30-m resolution continuous fields of tree cover: Landsat-based rescaling of MODIS vegetation continuous fields with lidar-based estimates of error. International Journal of Digital Earth, 6(5), 427-448.