In Nowosad/lopata: Local Pattern Analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  dev = "png"
)

The pattern-based spatial analysis makes it possible to search for areas with similar spatial patterns. This vignette shows how to do spatial patterns' search on example datasets. Let's start by attaching necessary packages:

library(motif)
library(stars)
library(sf)
library(tmap)

Spatial patterns' search requires two spatial objects. The first one is the area of interest, and the second one is a larger area that we want to search in. For this vignette, we read the "raster/landcover2015.tif" file, and crop our area of interest using coordinates of its borders.

landcover = read_stars(system.file("raster/landcover2015.tif", package = "motif"))

ext = st_bbox(c(xmin = 238000, xmax = 268000,
                ymin = -819814, ymax = -789814),
                crs = st_crs(landcover))

landcover_ext = landcover[ext]

The landcover_ext represents area mostly covered by forest with some agriculture.

landcover_ext = droplevels(landcover_ext)
plot(landcover_ext, key.pos = 4, key.width = lcm(5), main = NULL)

We want to compare it to the land cover dataset of New Guinea - landcover.

my_breaks = c(0, 0.01, 0.05, 0.1, 0.4, 0.7, 1.01)
my_palette = hcl.colors(n = 6, palette = "Vik")
landcover = droplevels(landcover)
plot(landcover, key.pos = 4, key.width = lcm(5), main = NULL)

Regular local landscapes

Spatial patterns' search is done by the lsp_search() function. It expects an area of interest as the first object and the larger area as the second one. We should provide the type of signature (type) and the suitable distance function (dist_fun) that we want to use to compare two datasets. Additional arguments include the size of the search window from the larger area (window) and how much of NA values we can accept in the local landscapes (threshold).

search_1 = lsp_search(landcover_ext, landcover, 
                      type = "cove", dist_fun = "jensen-shannon",
                      window = 100, threshold = 1)
search_1

The result of the lsp_search() function is a stars object with three attributes:

• id - unique value for each local landscape
• na_prop - proportion (0-1) of NA values in each local landscapes
• dist - the distance between our area of interest and each local landscape

We can visualize the results, using, for example, the tmap package:

my_breaks = c(0, 0.001, 0.01, 0.1, 1.01)
tm_shape(search_1) +
  tm_raster("dist", breaks = my_breaks, palette = "-viridis") +
  tm_layout(legend.outside = TRUE)

It is now possible to see that there are several areas with a distance below 0.001 represented by a yellow color - they are the most similar to landcover_ext.

min_dist1 = min(search_1$dist, na.rm = TRUE)
min_dist1

We can find their ids using the code below.

unique(search_1$id[which(search_1$dist < 0.001)])

To extract selected local landscape, the lsp_extract() function can be used.

search_1_690 = lsp_extract(landcover,
                            window = 100,
                            id = 690)

Its output is a stars object, that we can vizualize and see that it is fairly similar to the area of interest.

search_1_690 = droplevels(search_1_690)
plot(search_1_690, main = NULL)

Irregular local landscapes

Search is also possible in irregular local landscapes, based on polygon data. ecoregions.gpkg contains terrestrial ecoregions for New Guinea from https://ecoregions2017.appspot.com/.

ecoregions = read_sf(system.file("vector/ecoregions.gpkg", package = "motif"))

This dataset has 22 rows, where each row relates to one ecoregion. Each ecoregion is also related to a unique value in the id column.

ecoregions = st_transform(ecoregions, st_crs(landcover))

# https://medialab.github.io/iwanthue/
my_pal = c("#d34359", "#62b93c", "#b75fcf", "#53c069", "#d44295", "#acb939", "#626edd", "#dc9e36", "#8156a8", "#4b8734", "#d98dc7", "#58c096", "#cc542a", "#48bbd2", "#bf814d", "#6686c8", "#968c30", "#a34d78", "#36815b", "#c26963", "#a2b36b", "#6b6829")
plot(ecoregions["id"], main = NULL, col = my_pal, key.pos = 4, key.width = lcm(5))

The lsp_search() function works very similarly to the previous case - we just need to provide our ecoregions in the window argument.

search_2 = lsp_search(landcover_ext, landcover, 
                      type = "cove", dist_fun = "jensen-shannon",
                      window = ecoregions["id"], threshold = 1)
search_2

Let's vizualize the output:

my_breaks = c(0, 0.001, 0.01, 0.1, 1.01)
tm_shape(search_2) +
  tm_raster("dist", breaks = my_breaks, palette = "-viridis") +
  tm_shape(ecoregions) +
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE)

This search shows that most of the polygons are fairly different from our area of interest. Only one of them, located in the east, has a relatively small distance of about 0.007.

min_search2 = min(search_2$dist, na.rm = TRUE)
min_search2

We can obtain its id (10) using the code below.

unique(search_2$id[which(search_2$dist == min_search2)])

Now, we can use lsp_extract() to extract land cover for this ecoregion.

search_2_10 = lsp_extract(landcover,
                          window = ecoregions["id"],
                          id = 10)

This local landscape is also mostly covered by forest with just some smaller areas of agriculture.

search_2_10 = droplevels(search_2_10)
plot(search_2_10, main = NULL)

Nowosad/lopata documentation built on April 15, 2024, 4:32 p.m.