starsExtra: Miscellaneous Functions for Working with 'stars' Rasters

Description

Miscellaneous functions for working with 'stars' objects, mainly single-band rasters. Currently includes functions for: (1) focal filtering, (2) detrending of Digital Elevation Models, (3) calculating flow length, (4) calculating the Convergence Index, (5) calculating topographic aspect and topographic slope.

Author(s)

Maintainer: Michael Dorman dorman@post.bgu.ac.il

See Also

Useful links:

• https://michaeldorman.github.io/starsExtra/

• https://github.com/michaeldorman/starsExtra/

• Report bugs at https://github.com/michaeldorman/starsExtra/issues/

Calculate the Convergence Index (CI) from a slope raster

Description

Calculates the Convergence Index (CI) given a topographic slope raster. Input and output are rasters of class stars, single-band (i.e., only '"x"' and '"y"' dimensions), with one attribute.

Usage

CI(x, k, na.rm = FALSE, na_flag = -9999)

Arguments

Value

A stars raster with CI values.

Note

The raster is "padded" with (k-1)/2 more rows and columns of NA values on all sides, so that the neighborhood of the outermost rows and columns is still a complete neighborhood. Those rows and columns are removed from the final result before returning it. Aspect values of -1, specifying flat terrain, are assigned with a CI value of 0 regardless of their neighboring values.

References

The Convergence Index algorithm is described in:

Thommeret, N., Bailly, J. S., & Puech, C. (2010). Extraction of thalweg networks from DTMs: application to badlands.

Examples

# Small example
data(dem)
dem_asp = aspect(dem)
dem_ci = CI(dem_asp, k = 3)
r = c(dem, round(dem_ci, 1), along = 3)
r = st_set_dimensions(r, 3, values = c("input (aspect)", "output (CI, k=3)"))
plot(r, text_values = TRUE, breaks = "equal", col = terrain.colors(10), mfrow = c(1, 2))

# Larger example
data(golan)
golan_asp = aspect(golan)
golan_ci = CI(golan_asp, k = 25)
plot(golan_asp, breaks = "equal", col = hcl.colors(11, "Spectral"), main = "input (aspect)")
plot(golan_ci, breaks = "equal", col = hcl.colors(11, "Spectral"), main = "output (CI, k=25)")

Calculate topographic aspect from a DEM

Description

Calculates topographic aspect given a Digital Elevation Model (DEM) raster. Input and output are rasters of class stars, single-band (i.e., only '"x"' and '"y"' dimensions), with one attribute.

Usage

aspect(x, na_flag = -9999)

Arguments

Value

A stars raster with topographic slope, i.e., the azimuth where the terrain is tilted towards, in decimal degrees (0-360) clockwise from north. Aspect of flat terrain, i.e., where all values in the neighborhood are equal, is set to -1. Returned raster values are of class units (decimal degrees).

Note

Aspect calculation results in NA when at least one of the cell neighbors is NA, including the outermost rows and columns. Given that the focal window size in aspect calculation is 3*3, this means that the outermost one row and one column are given an aspect value of NA.

The raster must be in projected CRS, and units of x/y resolution are assumed to be the same as units of elevation (typically meters).

References

The topographic aspect algorithm is based on the How aspect works article in the ArcGIS documentation:

https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-aspect-works.htm

Examples

# Small example
data(dem)
dem_aspect = aspect(dem)
plot(
  dem, text_values = TRUE, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "input (elevation)"
)
plot(
  dem_aspect, text_values = TRUE, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "output (aspect)"
)

# Larger example
data(carmel)
carmel_aspect = aspect(carmel)
plot(
  carmel, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "input (elevation)"
)
plot(
  carmel_aspect, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "output (aspect)"
)

Digital Elevation Model of Mount Carmel

Description

A stars object representing a Digital Elevation Model (DEM) Digital Elevation Model of Mount Carmel, at 90m resolution

Usage

carmel

Format

A stars object with 1 attribute:

elevation

Elevation above sea level, in meters

Examples

plot(carmel, breaks = "equal", col = terrain.colors(11))

Small Digital Elevation Model

Description

A stars object representing a small 13*11 Digital Elevation Model (DEM), at 90m resolution

Usage

dem

Format

A stars object with 1 attribute:

elevation

Elevation above sea level, in meters

Examples

plot(dem, text_values = TRUE, breaks = "equal", col = terrain.colors(11))

Detrend a Digital Elevation Model

Description

Detrends a Digital Elevation Model (DEM) raster, by subtracting a trend surface. The trend is computed using mgcv::gam or mgcv::bam (when parallel>1) with formula z ~ s(x, y).

Usage

detrend(x, parallel = 1)

Arguments

Value

A two-dimensional stars object, with two attributes:

• resid - the detrended result, i.e., "residual"

• trend - the estimated "trend" which was subtracted from the actual elevation to obtain resid

Examples

# Small example
data(dem)
dem1 = detrend(dem)
dem1 = st_redimension(dem1)
dem1 = st_set_dimensions(dem1, 3, values = c("resid", "trend"))
plot(round(dem1), text_values = TRUE, col = terrain.colors(11))

# Larger example 1
data(carmel)
carmel1 = detrend(carmel, parallel = 2)
carmel1 = st_redimension(carmel1)
carmel1 = st_set_dimensions(carmel1, 3, values = c("resid", "trend"))
plot(carmel1, col = terrain.colors(11))

# Larger example 2
data(golan)
golan1 = detrend(golan, parallel = 2)
golan1 = st_redimension(golan1)
golan1 = st_set_dimensions(golan1, 3, values = c("resid", "trend"))
plot(golan1, col = terrain.colors(11))

Calculate raster of distances to nearest feature

Description

Given a stars raster and an sf vector layer, returns a new raster with the distances of each cell centroid to the nearest feature in the vector layer.

Usage

dist_to_nearest(x, v, progress = TRUE)

Arguments

Value

A stars raster with distances to nearest feature

Examples

# Sample 'sf' layer
x = st_point(c(0,0))
y = st_point(c(1,1))
x = st_sfc(x, y)
x = st_sf(x)
x = st_buffer(x, 0.5)

# Make grid
r = make_grid(x, res = 0.1, buffer = 0.5)
d = dist_to_nearest(r, x, progress = FALSE)

# Plot
plot(d, breaks = "equal", axes = TRUE, reset = FALSE)
plot(st_geometry(x), add = TRUE, pch = 4, cex = 3)

Extract raster values by lines or polygons

Description

Extract raster values by lines or polygons, summarizing for each feature using a function specified by the user. This function is aimed to reproduce (some of) the functionality of raster::extract.

Usage

extract2(x, v, fun, progress = TRUE, ...)

Arguments

Value

A vector (single-band raster) or matrix (multi-band raster) with the extracted and summarized values

Examples

# Polygons
pol = st_bbox(landsat)
pol = st_as_sfc(pol)
set.seed(1)
pol = st_sample(pol, 5)
pol = st_buffer(pol, 100)
pol = c(pol, pol)

# Plot
plot(landsat[,,,1,drop=TRUE], reset = FALSE)
plot(pol, add = TRUE)

# Single-band raster
aggregate(landsat[,,,1,drop=TRUE], pol, mean, na.rm = TRUE)[[1]]  ## Duplicated areas get 'NA'
extract2(landsat[,,,1,drop=TRUE], pol, mean, na.rm = TRUE, progress = FALSE)

# Multi-band example
extract2(landsat, pol, mean, na.rm = TRUE, progress = FALSE)

# Lines
lines = st_cast(pol, "LINESTRING")

# Single-band raster
extract2(landsat[,,,1,drop=TRUE], lines, mean, na.rm = TRUE, progress = FALSE)

# Multi-band example
extract2(landsat, lines, mean, na.rm = TRUE, progress = FALSE)

Calculate flow length

Description

Calculates flow length for each pixel in a Digital Elevation Model (DEM) raster. Inputs and output are rasters of class stars, single-band (i.e., only '"x"' and '"y"' dimensions), with one attribute.

Usage

flowlength(elev, veg, progress = TRUE)

Arguments

Value

A numeric stars raster where each cell value is flow length, in resolution units.

References

The algorithm is described in:

Mayor, A. G., Bautista, S., Small, E. E., Dixon, M., & Bellot, J. (2008). Measurement of the connectivity of runoff source areas as determined by vegetation pattern and topography: A tool for assessing potential water and soil losses in drylands. Water Resources Research, 44(10).

Examples

# Example from Fig. 2 in Mayor et al. 2008

elev = rbind(
  c(8, 8, 8, 8, 9, 8, 9),
  c(7, 7, 7, 7, 9, 7, 7),
  c(6, 6, 6, 6, 6, 5, 7),
  c(4, 5, 5, 3, 5, 4, 7),
  c(4, 5, 4, 5, 4, 6, 5),
  c(3, 3, 3, 3, 2, 3, 3),
  c(2, 2, 2, 3, 4, 1, 3)
)
veg = rbind(
  c(TRUE,  TRUE,  TRUE,  TRUE,  FALSE, FALSE, TRUE),
  c(TRUE,  TRUE,  TRUE,  TRUE,  TRUE,  FALSE, FALSE),
  c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE),
  c(FALSE, TRUE,  FALSE, FALSE, FALSE, FALSE, TRUE),
  c(TRUE,  TRUE,  FALSE, FALSE, FALSE, FALSE, FALSE),
  c(TRUE,  TRUE,  TRUE,  FALSE, FALSE, FALSE, FALSE),
  c(FALSE, TRUE,  TRUE,  FALSE, FALSE, TRUE,  TRUE)
)
elev = matrix_to_stars(elev)
veg = matrix_to_stars(veg)

# Calculate flow length
fl = flowlength(elev, veg, progress = FALSE)

# Plot
plot(
  round(elev, 1), text_values = TRUE, breaks = "equal", 
  col = terrain.colors(6), main = "input (elevation)"
)
plot(
  veg*1, text_values = TRUE, breaks = "equal", 
  col = rev(terrain.colors(2)), main = "input (vegetation)"
)
plot(
  round(fl, 1), text_values = TRUE, breaks = "equal", 
  col = terrain.colors(6), main = "output (flowlength)"
)

# Larger example
data(carmel)
elev = carmel
elev[is.na(elev)] = 0
veg = elev > 100
fl = flowlength(elev, veg, progress = FALSE)
plot(fl, breaks = "equal", col = hcl.colors(11), main = "flowlength (m)")

Apply a focal filter on a raster

Description

Applies a focal filter with weighted neighborhood w on a raster. The weights (w) can be added to, subtracted from, multiplied by or divided with the raster values (as specified with weight_fun). The focal cell is then taken as the mean, sum, minimum or maximum of the weighted values (as specified with fun). Input and output are rasters of class stars, single-band (i.e., only '"x"' and '"y"' dimensions), with one attribute.

Usage

focal2(
  x,
  w,
  fun = "mean",
  weight_fun = "*",
  na.rm = FALSE,
  mask = FALSE,
  na_flag = -9999
)

Arguments

Value

The filtered stars raster.

Note

The raster is "padded" with (nrow(w)-1)/2 more rows and columns of NA values on all sides, so that the neighborhood of the outermost rows and columns is still a complete neighborhood. Those rows and columns are removed from the final result before returning it. This means, for instance, that the outermost rows and columns in the result will be NA when using na.rm=FALSE.

References

The function interface was inspired by function raster::focal. The C code for this function is a modified and expanded version of the C function named applyKernel included with R package spatialfil.

Examples

# Small example
data(dem)
dem_mean3 = focal2(dem, matrix(1, 3, 3), "mean")
r = c(dem, round(dem_mean3, 1), along = 3)
r = st_set_dimensions(r, 3, values = c("input", "output (mean, k=3)"))
plot(r, text_values = TRUE, breaks = "equal", col = terrain.colors(11))

# Larger example
data(carmel)
carmel_mean15 = focal2(carmel, matrix(1, 15, 15), "mean")
r = c(carmel, carmel_mean15, along = 3)
r = st_set_dimensions(r, 3, values = c("input", "output (mean, k=15)"))
plot(r, breaks = "equal", col = terrain.colors(11))

Apply a focal filter on a raster (R)

Description

Applies a focal filter with neighborhood size k*k on a raster (class stars), using R code. This function is relatively slow, provided here mainly for testing purposes or for custom using functions which are not provided by focal2.

Usage

focal2r(x, w, fun, mask = FALSE, ...)

Arguments

Value

The filtered stars raster

Note

The raster is "padded" with one more row/column of NA values on all sides, so that the neigborhood of the outermost rows and columns is still a complete 3x3 neighborhood. Those rows and columns are removed from the filtered result before returning it.

Examples

# Small example
data(dem)
dem1 = focal2r(dem, matrix(1,3,3), mean, na.rm = TRUE)
dem2 = focal2r(dem, matrix(1,3,3), min, na.rm = TRUE)
dem3 = focal2r(dem, matrix(1,3,3), max, na.rm = TRUE)
r = c(dem, round(dem1, 1), dem2, dem3, along = 3)
r = st_set_dimensions(r, 3, values = c("input", "mean", "min", "max"))
plot(r, text_values = TRUE, breaks = "equal", col = terrain.colors(10))
# Larger example
data(carmel)
carmel1 = focal2r(carmel, matrix(1,3,3), mean, na.rm = TRUE, mask = TRUE)
carmel2 = focal2r(carmel, matrix(1,9,9), mean, na.rm = TRUE, mask = TRUE)
carmel3 = focal2r(carmel, matrix(1,15,15), mean, na.rm = TRUE, mask = TRUE)
r = c(carmel, carmel1, carmel2, carmel3, along = 3)
r = st_set_dimensions(r, 3, values = c("input", "k=3", "k=9", "k=15"))
plot(r, breaks = "equal", col = terrain.colors(100))

Footprints

Description

Calculates a polygon layer with the footprints of raster images.

Usage

footprints(x)

Arguments

Value

An sf layer with the footprints (i.e., bounding box polygons) of the rasters

Examples

# Create sample files
file1 = tempfile(fileext = ".tif")
file2 = tempfile(fileext = ".tif")
file3 = tempfile(fileext = ".tif")
r1 = landsat[,1:100, 1:100,]
r2 = landsat[,101:200, 101:200,]
r3 = landsat[,21:40, 51:120,]
write_stars(r1, file1)
write_stars(r2, file2)
write_stars(r3, file3)

# Calculate footprints
files = c(file1, file2, file3)
pol = footprints(files)
pol

Digital Elevation Model of Golan Heights

Description

A stars object representing a Digital Elevation Model (DEM) Digital Elevation Model of part of the Golan Heights and Lake Kinneret, at 90m resolution

Usage

golan

Format

A stars object with 1 attribute:

elevation

Elevation above sea level, in meters

Examples

plot(golan, breaks = "equal", col = terrain.colors(11))

RGB image of Mount Carmel

Description

A stars object representing an RGB image of part of Mount Carmel, at 30m resolution. The data source is Landsat-8 Surface Reflectance product.

Usage

landsat

Format

A stars object with 1 attribute:

refl

Reflectance, numeric value between 0 and 1

Examples

plot(landsat, breaks = "equal")

Get stars layer values as matrix

Description

Extracts the values of a single layer in a stars object to a matrix.

Usage

layer_to_matrix(x, check = TRUE)

Arguments

Value

A matrix with the layer values, having the same orientation as the raster (i.e., rows represent the y-axis and columns represent the x-axis).

Examples

data(dem)
m = layer_to_matrix(dem)
m

Get stars layer values as vector

Description

Extracts the values of a single layer in a stars object to a vector. Cell values are ordered from top-left corner to the right.

Usage

layer_to_vector(x, check = TRUE)

Arguments

Value

A vector with cell values, ordered by rows, starting from the top left corner (north-west) and to the right.

Examples

data(dem)
v = layer_to_vector(dem)
v

Make 'stars' grid from 'sf' layer

Description

Create 'stars' raster grid from bounding box of 'sf' vector layer, possibly buffered, with specified resolution.

Usage

make_grid(x, res, buffer = 0)

Arguments

Value

A stars raster with the grid, with all cell values equal to 1

Examples

# Sample 'sf' layer
x = st_point(c(0,0))
y = st_point(c(1,1))
x = st_sfc(x, y)
x = st_sf(x)

# Make grid
r = make_grid(x, res = 0.1, buffer = 0.5)
r[[1]][] = rep(1:3, length.out = length(r[[1]]))

# Plot
plot(r, axes = TRUE, reset = FALSE)
plot(st_geometry(x), add = TRUE, pch = 4, cex = 3, col = "red")

Extend matrix

Description

Adds n rows and columns with NA values on all sides of a matrix.

Usage

matrix_extend(m, n = 1, fill = NA)

Arguments

Value

An extended matrix

Examples

m = matrix(1:6, nrow = 2, ncol = 3)
m
matrix_extend(m, 1)
matrix_extend(m, 2)
matrix_extend(m, 3)

Get neighboring cell values for given matrix cell

Description

Get the values of a k*k neighborhood, as vector and by row, given a matrix, k, and focal cell position (row and column).

Usage

matrix_get_neighbors(m, pos, k = 3)

Arguments

Value

A vector with cell values, ordered by rows, starting from the top left corner of the neighborhood and to the right. When neighborhood extends beyond matrix bounds, only the "existing" values are returned.

Examples

m = matrix(1:12, nrow = 3, ncol = 4)
m
matrix_get_neighbors(m, pos = c(2, 2), k = 3)
matrix_get_neighbors(m, pos = c(2, 2), k = 5)
matrix_get_neighbors(m, pos = c(2, 2), k = 7)  # Same result

Convert matrix to stars

Description

Converts matrix to a single-band stars raster, conserving the matrix orientation where rows become the y-axis and columns become the y-axis. The bottom-left corner of the axis is set to (0,0) coordinate, so that x and y coordinates are positive across the raster extent.

Usage

matrix_to_stars(m, res = 1)

Arguments

Value

A stars raster

Examples

data(volcano)
r = matrix_to_stars(volcano, res = 10)
plot(r)

Trim matrix

Description

Removes n rows and columns with NA values on all sides of a matrix.

Usage

matrix_trim(m, n = 1)

Arguments

Value

A trimmed matrix, or NULL if trimming results in an "empty" matrix.

Examples

m = matrix(1:80, nrow = 8, ncol = 10)
m
matrix_trim(m, 1)
matrix_trim(m, 2)
matrix_trim(m, 3)
matrix_trim(m, 4)

Mode

Description

Find the mode (i.e., most common value) in a numeric vector. In case of ties, the first encountered value is returned.

Usage

mode_value(x, na_flag = -9999)

Arguments

Value

The mode, numeric vector of length 1

Examples

x = c(3, 2, 5, 5, 3, 10, 2, 5)
mode_value(x)

Normalize a 2D 'stars' object

Description

Check, and possibly correct, that the input stars object:

• Has exactly one attribute.

• Has exactly two dimensions.

• The dimensions are spatial, named x and y (in that order).

Usage

normalize_2d(x)

Arguments

Value

A new stars object

Examples

# Small example
data(dem)
normalize_2d(dem)

Normalize a 3D 'stars' object

Description

Check, and possibly correct, that the input stars object:

• Has exactly one attribute.

• Has exactly three dimensions.

• The first two dimensions are spatial, named x and y (in that order).

Usage

normalize_3d(x)

Arguments

Value

A new stars object

Examples

# Small example
data(landsat)
normalize_3d(landsat)

Convert RGB to greyscale

Description

Convert a 3-band RGB raster to 1-band greyscale raster. Based on wvtool::rgb2gray.

Usage

rgb_to_greyscale(x, rgb = 1:3, coefs = c(0.3, 0.59, 0.11))

Arguments

Value

A two-dimensional stars object greyscale values

Examples

data(landsat)
plot(landsat, rgb = 1:3)
landsat_grey = rgb_to_greyscale(landsat)
plot(landsat_grey, breaks = "equal")

Calculate topographic slope from a DEM

Description

Calculates topographic slope given a Digital Elevation Model (DEM) raster. Input and output are rasters of class stars, single-band (i.e., only '"x"' and '"y"' dimensions), with one attribute.

Usage

slope(x, na_flag = -9999)

Arguments

Value

A stars raster with topographic slope, i.e., the azimuth where the terrain is tilted towards, in decimal degrees (0-360) clockwise from north.

Note

Slope calculation results in NA when at least one of the cell neighbors is NA, including the outermost rows and columns. Given that the focal window size in slope calculation is 3*3, this means that the outermost one row and one column are given an slope value of NA.

The raster must be in projected CRS, and units of x/y resolution are assumed to be the same as units of elevation (typically meters).

References

The topographic slope algorithm is based on the How slope works article in the ArcGIS documentation:

https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-slope-works.htm

Examples

# Small example
data(dem)
dem_slope = slope(dem)
plot(
  dem, text_values = TRUE, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "input (elevation)"
)

plot(
  dem_slope, text_values = TRUE, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "output (slope)"
)
# Larger example
data(carmel)
carmel_slope = slope(carmel)
plot(
  carmel, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "input (elevation)"
)
plot(
  carmel_slope, breaks = "equal", 
  col = hcl.colors(11, "Spectral"), main = "output (slope)"
)

Remove empty outer rows and columns

Description

Removes complete outer rows and columns which have NA values.

Usage

trim2(x)

Arguments

Value

A new stars object with empty outer rows and columns removed

Examples

# Single-band example
data(dem)
dem[[1]][1,] = NA
dem1 = trim2(dem)

# Multi-band example
data(landsat)
landsat[[1]][1:100,,] = NA
landsat1 = trim2(landsat)

Create matrix with azimuths to center

Description

Creates a matrix with directions (i.e., azimuth) to central cell, of specified size k. The matrix can be used as weight matrix when calculating the convergence index (see Examples).

Usage

w_azimuth(k)

Arguments

Value

A matrix where each cell value is the azimuth from that cell towards the matrix center.

Examples

m = w_azimuth(3)
m
m = w_azimuth(5)
m

Create matrix with circular weight pattern

Description

Creates a matrix with where a circular pattern is filled with values of 1 and the remaining cells are filled with values of 0 (see Examples).

Usage

w_circle(k)

Arguments

Value

A matrix with a circular pattern.

Examples

m = w_circle(3)
image(m, asp = 1, axes = FALSE)
m = w_circle(5)
image(m, asp = 1, axes = FALSE)
m = w_circle(15)
image(m, asp = 1, axes = FALSE)
m = w_circle(35)
image(m, asp = 1, axes = FALSE)
m = w_circle(91)
image(m, asp = 1, axes = FALSE)
m = w_circle(151)
image(m, asp = 1, axes = FALSE)