Getting started with landscapeDecoupler

A little about the concept

The image bellow summarizes the main ideas about coupled/decoupled multiscale approaches. You may noticed that the “traditional” nested landscape approach always have the other scales inside it. So, one may be interested in see how composition and structure of each portion of the landscape can contribute to any variable responses. Therefore, users can either use a sequential decoupling landscape strategy (Figure 1, B) or cut out an specific scale and compare against the total (Figure 1, C and D).

Figure 1. A summary of multiscale strategies

The idea of the LandscapeDecoupler is to recursively extract and analyze decoupled buffered scales from a given landscape, returning raster objects, files and basic metrics of the landscape components, so users can easily use rasters and metrics in further analysis.

The package needs three inputs:
1. a classified landscape (a raster)
2. lat-long coordinates (a txt or shape [.shp] file);
3. two or more buffer sizes to be decoupled (input as vector).

The package also includes some conveniences to plot and export the results.

Sequential decoupling

We first load the package:

library(landscapeDecoupler)
library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(raster)

## Loading required package: sp

Then the landscape:
You probably already have a classified raster, if so, you just need to raster() from its path…
Herein we will load an example from MapBiomas 5.0 collections.

#1 - loading some raster
file <- system.file("extdata/raster.grd", package="landscapeDecoupler")
r <- raster(file)
#1.1 - making categories
r <- as.factor(r)

Grabbing some meaningful colors
The package provide several color pallets from the grDevices, but we will load an specific index of colors used by MapBiomas project. This file is a table with corresponding the raster values and colors.

#2 - loading the color index file
color.file <- system.file("extdata/colors.txt",
                    package="landscapeDecoupler")
colors <- read.table( color.file , header = T, sep = "\t")
#How it looks like
head(colors)

##   value              name    R    G    B Alpha     hex
## 1     1            Forest 0.07 0.60 0.07     1 #129912
## 2     2       Nat. Forest 0.12 0.27 0.14     1 #1F4524
## 3     3      Forest Form. 0.00 0.39 0.00     1 #006300
## 4     4     Savanna Form. 0.00 1.00 0.00     1 #00FF00
## 5     5          Mangrove 0.41 0.46 0.22     1 #697538
## 6     9 Forest Plantation 0.58 0.32 0.20     1 #945233

We now find the unique values from the raster and use it to filter only the colors presented in our landscape.

#2.1 - creating some breaks to match categories and colors
breaks <- c(as.numeric(levels(as.factor(getValues(r)))))
#2.2 - matching index values and breaks
colors <- colors[colors$value %in% breaks,]
#2.3 - ordering by values
colors <- colors[order(colors$value),]

We can check it by plotting.

# A simple plot to see the landscape and categories meaning.
par(mar=c(3,2,3,11))
plot(r, col= colors$hex, legend=F, 
     xpd=T, main="An eight classes raster")
legend(x=250000, y=7725000, 
       legend = paste0(colors$value, ": ", colors$name),
       title = "Values / Categories:", title.adj = 0.05,
       fill = colors$hex, bty = "n", xpd = T)

Now to the points

Users can read latitude and longitude files through the sf package functions. Users are encouraged to provide an “id” field in shape files in order to keep custom site names. Otherwise, sites will be sequentially named.

In this example we will load an shp file with the function read_points.

IMPORTANT shape files or txt files must have a field named “id” as it will be used as unique identifier for the outputs

#2 - reading points
points.file <- system.file("extdata/pnts.shp", package="landscapeDecoupler")
p <- read_points(points.file, type = "shp")
plot(r, col=colors$hex, legend=FALSE)
plot(p, pch=19, col="black", add=T)
text(x=p, p$id, adj=c(0.5, -0.5), font=2, col="black", cex=1.1)

In any analysis of spatial data it’s a good practice to keep coordinate systems consistent. Shape files usually already have a crs assigned, but users might be shure using the code bellow.

crs(p)<-crs(r)

And buffers:

Users must provide a set of buffers (> 2) to be cropped out of the landscape. It must be provided in a sequential and increasing order, as the example bellow.

#Setting the buffer sizes
b <- c(250, 500, 1000, 1500)

#To plot buffers
ps  <- sf::st_as_sf(do.call(rbind, replicate(4,p)))
bfs <- sf::st_buffer(ps, rep(b, 4))
#calling the plots
plot(r, col=colors$hex)
points(p, pch=16, col="blue", cex=0.5)
plot(bfs$geometry, add=T)

Users can provide any set of buffers, but you must bear in mind that the final landscapes (scales) are calculated as the difference between sequential buffers (excepting the one corresponding to the core, usually the first). Thus, the first scale (250) will be a circle with 250 m radius from the sites we provided, but the second scale (500) will be the 500 m radius circle less the first buffer (250); Likewise, the third buffer, 1000 m less the 500 m, and so on through the list.

Now we can run the main function of the package, asking for decoupling a landscape r at a site p with a set of buffers b.

multiscales <- decouple(r, p, b)

## Loading required namespace: rgeos

It might take a couple seconds or minuts depending on the resolution of the raster, the number of sites and buffers.

We can take a look at the object returned. It is a list of several RasterBrick objects ( p x b). Each item on the list correspond to a sampling site and each scale is a layer in the RasterBrick

Rasters on the list are named following the site ID, as provided by the p object (.shp file).

It’s very streightforward to access any scale and sampling site. The sintax multiscales$p01$X250 will return the 250 metters buffer from the p01 site.

multiscales$p01$X250

## class      : RasterLayer 
## dimensions : 99, 108, 10692  (nrow, ncol, ncell)
## resolution : 27.59161, 30.31018  (x, y)
## extent     : 198126.4, 201106.3, 7680610, 7683611  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=24 +south +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : X250 
## values     : 3, 3  (min, max)
## attributes :
##        ID
##  from:  3
##   to : 41

traditional list sintax will also works. The following sintax will return the same: multiscales[[1]][[1]]

We can extract very basic information of the raster with the function extract_metrics just like bellow:

metrics <- extract_metrics(multiscales, countNA = "no")
head(metrics)

##   site layer class metric value
## 1  p01  X250     3 counts   239
## 2  p01  X500     3 counts   667
## 3  p01  X500    15 counts    15
## 4  p01  X500    21 counts    17
## 5  p01 X1000     3 counts  1445
## 6  p01 X1000    15 counts   807

The output from extract_metrics is a data frame with very basic stats, as class counting, its percentage, total number of classes in a given scale and compositional heterogeneity.

Finally we can plot it. The package has a wrapper for the native plot function which will allow us to plot some sets of scales directly.

#to plot the first three sites
multiplot(multiscales[1:2])

## [[1]]

## 
## [[2]]

To plot scales from a particular site, changing the color pallet…

multiplot(multiscales$p01, palette = "viridis")

Plotting with custom colors.

#p01 and p02 plots in a 4x1 grid
multiplot(multiscales[1:2], colors = colors$hex,
          extra.arg = list(layout=c(4,1)))

## [[1]]

## 
## [[2]]

At last, we can save the decoupled scales in local files to access in other software or packages.

write_scales(x = multiscales, path = getwd(), format = 'GTiff')

Specific decoupling

Users might want to decouple specific scales from a “doucoupled” landscape, for instance, decouple a scale between 250 and 500 m out of a landscape with 1500 m radius total. In the LandscapeDecouple package it is simply done as bellow:

specific.scales <- decouple.specific(multiscales, specific = 2)
#you can also decouple more scales at once, like bellow:
#specific.scales <- decouple.specific(multiscales, c(2,4))

for visual inspection, the user can plot the resulted scales at the site 01:

multiplot(specific.scales[[1]], colors = colors$hex)

It is easy to see that we are fully able to compare any information from one scale against to the whole left. You may also notice that the name of the “nonspecific” scale is a concatenation of two other buffers we used to configure our landscape.

If we extract the metrics from these two raster it will seems like this:

specific.scales.metrics <- extract_metrics(specific.scales)

head(specific.scales.metrics)

##   site         layer class metric value
## 1  p01          X500     3 counts   667
## 2  p01          X500    15 counts    15
## 3  p01          X500    21 counts    17
## 4  p01 merged_scales     3 counts  3528
## 5  p01 merged_scales     9 counts    12
## 6  p01 merged_scales    15 counts  2409

extract_metrics is a very fast and convinient way to get some metrics from two different levels: classes and landscape counts = number of cells for a given class heterog = the landscape heterogeneity n_class = the number of classes in a given landscape percent = the proportion of each class in a given landscape

The package also has the possibility to access any metrics in landscapeMetrics package. See details in following sections.

Ploting with other packages

Users can plot decouple output with different packages, including R base and tmap, for instance…

You just need to stack the list or make sure you are pointing to an specific raster item on the list.

For R base / raster package plot

#For R base
plot(specific.scales[[1]][[2]])

#or plot(specific.scales$p01$X500)
#or plot(multiscales[[1]][[4]])
#or plot(multiscales$p01$X250)

IMPORTANT if you “id” field is not alphanumeric R will fale in use the $ operator to access list names. Therefore, you will have to use an syntax like myobject[1] to get accesses to your scales in the same order you passed it to the function.

Plotting with tmap

tmap is a quite popular map plotting package and can be easily used to produce more elaborated plots:

library("tmap")

## Registered S3 methods overwritten by 'stars':
##   method             from
##   st_bbox.SpatRaster sf  
##   st_crs.SpatRaster  sf

#creating a map for the site p01
map1 <- specific.scales$p01
#checking the categories in this stack
mapcat <- levels(as.factor(getValues(map1)))
#filtering colors that match with stacks
mapcols <- colors[ colors$value %in% mapcat ,]
#now plotting
tm_shape(map1) +
  tm_raster(style = "cat",
            pal = mapcols$hex,
            title = "Values / Categories:",
            labels = paste0(mapcols$value, ": ", mapcols$name)) +
  tm_compass() +
  tm_scale_bar() +
  tm_grid(n.y=3, n.x=3, lines = FALSE, labels.rot = c(0, 90)) +
  tm_layout(legend.outside = T, 
            legend.outside.position = "bottom")

Calculating Metrics

Once the user has defined the desired level of comparison and extract his landscape, one might want to proceed with more complex analysis. Other packages as LandscapeMetrics might be used to calculate landscape properties.

the landscapeDecoupler has an wrapper for the main LandscapeMetrics function, the calculate_lsm. You can call calculations from a desired metric implemented in landscapeMetrics as follow:

#creating an object to receive the data frame
metrics.from.lsm <- calc_lsm(multiscales, level="class", metric=c("ai","ca"))
#printing the head of the data frame
head(metrics.from.lsm)

## # A tibble: 6 × 7
##   site  layer level class id    metric value
##   <fct> <fct> <fct> <fct> <fct> <fct>  <dbl>
## 1 p01   X250  class 3     <NA>  ai     99.3 
## 2 p01   X250  class 3     <NA>  ca     20.0 
## 3 p01   X500  class 3     <NA>  ai     95.9 
## 4 p01   X500  class 3     <NA>  ca     55.8 
## 5 p01   X500  class 15    <NA>  ai     68.2 
## 6 p01   X500  class 15    <NA>  ca      1.25

As you can see, it is quite easy to access any metrics from landscapeMetrics using calc_lsm. You can list available metrics with list_lsm().

library("landscapemetrics")
area <- calculate_lsm(specific.scales$p01, level = "class", metric = "area")
#Just to show you who are layers 1 and 2 for the landscapeMetrics
paste0("layer", c(1:2),": ", names(specific.scales$p01))

## [1] "layer1: X500"          "layer2: merged_scales"

area[area$metric=="area_mn",  ]

## # A tibble: 7 × 6
##   layer level class    id metric   value
##   <int> <chr> <int> <int> <chr>    <dbl>
## 1     1 class     3    NA area_mn 55.8  
## 2     1 class    15    NA area_mn  1.25 
## 3     1 class    21    NA area_mn  0.711
## 4     2 class     3    NA area_mn 21.1  
## 5     2 class     9    NA area_mn  1.00 
## 6     2 class    15    NA area_mn 11.2  
## 7     2 class    21    NA area_mn  4.30

#You can also do like this:
#landscapes <- stack(specific.scales[[1]])
#area<-calculate_lsm(landscape, level = "class", metric = "area")
#area

lsmetrics_area <- lapply(specific.scales, calculate_lsm, "class", "area" )
head(lsmetrics_area$p01)

## # A tibble: 6 × 6
##   layer level class    id metric    value
##   <int> <chr> <int> <int> <chr>     <dbl>
## 1     1 class     3    NA area_cv  NA    
## 2     1 class    15    NA area_cv  NA    
## 3     1 class    21    NA area_cv 125.   
## 4     1 class     3    NA area_mn  55.8  
## 5     1 class    15    NA area_mn   1.25 
## 6     1 class    21    NA area_mn   0.711

#or head(lsmetrics_area$p02) #and so on...