| Title: | Methods for Calculating Gradient Surface Metrics |
|---|---|
| Description: | Methods for calculating gradient surface metrics for continuous analysis of landscape features. |
| Authors: | Annie C. Smith, Phoebe Zarnetske, Kyla Dahlin, Adam Wilson, Andrew Latimer |
| Maintainer: | Annie C. Smith <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.1.0 |
| Built: | 2025-03-12 05:15:38 UTC |
| Source: | https://github.com/bioxgeo/geodiv |
Calculates the various texture metrics over a window centered
on an individual pixel. This function is modified slightly from the
calculate_lsm_focal function in the landscapemetrics package (Hesselbarth et al. 2019).
.calculate_met_focal(landscape, n_row, n_col, points, what, ...).calculate_met_focal(landscape, n_row, n_col, points, what, ...)
landscape |
A raster or matrix. |
n_row |
Numeric. Number of rows in focal window. |
n_col |
Numeric. Number of columns in focal window. |
points |
Dataframe. Coordinates and values of cells, calculated with the *landscapemetrics*
|
what |
Character. Metric to calculate for each window. Metrics from the geodiv package are listed below. |
... |
Additional arguments for the metric functions. All applicable arguments will be applied to the entire list of metrics. |
Metrics from geodiv package:
'sa': average surface roughness
'sq': root mean square roughness
's10z': ten-point height
'sdq': root mean square slope of surface, 2-point method
'sdq6': root mean square slope of surface, 7-point method
'sdr': surface area ratio
'sbi': surface bearing index
'sci': core fluid retention index
'ssk': skewness
'sku': kurtosis
'sds': summit density
'sfd': 3d fractal dimension
'srw': dominant radial wavelength, radial wavelength index, mean half wavelength
'std': angle of dominating texture, texture direction index
'svi': valley fluid retention index
'stxr': texture aspect ratio
'ssc': mean summit curvature
'sv': maximum valley depth
'sph': maximum peak height
'sk': core roughness depth
'smean': mean peak height
'svk': reduced valley depth
'spk': reduced peak height
'scl': correlation length
'sdc': bearing area curve height interval
The metric value over the window.
Hesselbarth, M.H.K., Sciaini, M., With, K.A., Wiegand, K., Nowosad, J. 2019. landscapemetrics: an open-source R tool to calculate landscape metrics. - Ecography 42:1648-1657(ver. 0).
Converts degree value(s) to radians.
.deg2rad(x).deg2rad(x)
x |
Numeric. Degree value(s). |
Numeric of degree value(s).
Internal function to calculates the maximum distances to specified autocorrelation values (e.g., 0.2) of the areal autocorrelation function (AACF). All 180 degrees from the origin of the AACF image are considered for the calculation.
.maxdist(threshold, aacfimg, distimg).maxdist(threshold, aacfimg, distimg)
threshold |
A number with a value between 0 and 1. Indicates the autocorrelation value to which the rates of decline are measured. |
aacfimg |
A raster of the areal autocorrelation function. This is the AACF raster split in two in terms of height. |
distimg |
A raster of distances to all pixels from the center of the original image. Distances are in meters if original raster was unprojected, and are in map units (usually meters) if raster was projected (see raster::distance documentation for more details). |
A list containing the maximum distances from an autocorrelation value of 1 to the specified autocorrelation value < 1. Distances are meters if original raster was unprojected, and are in map units (usually meters) if raster was projected (see raster::distance documentation for more details).
Internal function to calculates the minimum distances to specified autocorrelation values (e.g., 0.2) of the areal autocorrelation function (AACF). All 180 degrees from the origin of the AACF image are considered for the calculation.
.mindist(threshold, aacfimg, distimg).mindist(threshold, aacfimg, distimg)
threshold |
A number with a value between 0 and 1. Indicates the autocorrelation value to which the rates of decline are measured. |
aacfimg |
A raster of the areal autocorrelation function. This is the AACF raster split in two in terms of height. |
distimg |
A raster of distances to all pixels from the center of the original image. Distances are in meters if original raster was unprojected, and are in map units (usually meters) if raster was projected (see raster::distance documentation for more details). |
A list containing the minimum distances from an autocorrelation value of 1 to the specified autocorrelation value < 1. Distances are meters if original raster was unprojected, and are in map units (usually meters) if raster was projected (see raster::distance documentation for more details).
Converts radian value(s) to degrees.
.rad2deg(x).rad2deg(x)
x |
Numeric. Radian value(s). |
Numeric of degree value(s).
Calculates the areal autocorrelation function (AACF) as the
inverse of the Fourier power spectrum. aacf(x) returns
the AACF in both matrix and raster format.
aacf(x)aacf(x)
x |
An n x n raster or matrix. |
A raster or matrix representation of the AACF. Both raster and matrix values are normalized so that the maximum is equal to 1.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate aacf img and matrix aacf_out <- aacf(normforest) # plot resulting aacf image terra::plot(aacf_out)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate aacf img and matrix aacf_out <- aacf(normforest) # plot resulting aacf image terra::plot(aacf_out)
Calculates the area above the bearing area curve from
points a to b. If a box is drawn around
a function with the upper-left at a and the
bottom-right at b, this function extracts the area
above the function within the box.
area_above(f, a, b, n = 100)area_above(f, a, b, n = 100)
f |
The function for the Bearing Area curve produced by
|
a |
Numeric. The left x boundary. |
b |
Numeric. The right x boundary. |
n |
Numeric. The number of subdivisions along the function line. |
The area under the curve used to calculate area above the
curve is calculated as the numerical integral
of the Bearing Area function from a to b
using the trapezoid rule with n subdivisions. Assume
a < b and n is a positive integer.
A numeric value representing the area above the curve with
x bounds a and b.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # basic values z <- terra::values(normforest) # calculate cumulative probability density function of surface 'height' (= ndvi) mod <- ecdf((1 - z)) # valley fluid retention index = void volume in 'valley' zone Svi <- area_above(f = mod, b = 1, a = 0.8, n = 500)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # basic values z <- terra::values(normforest) # calculate cumulative probability density function of surface 'height' (= ndvi) mod <- ecdf((1 - z)) # valley fluid retention index = void volume in 'valley' zone Svi <- area_above(f = mod, b = 1, a = 0.8, n = 500)
Finds a rotated version of the Bearing Area (Abbott-Firestone) curve from a raster or matrix. The resulting function should be rotated 90 degrees clockwise to get the actual Bearing Area curve.
bearing_area(x)bearing_area(x)
x |
A raster or matrix. |
A function describing the rotated Bearing Area curve.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the rotated Bearing Area curve. ba_func <- bearing_area(normforest) # rotate the values and re-plot xval <- environment(ba_func)$y yval <- (1 - environment(ba_func)$x) plot(yval ~ xval)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the rotated Bearing Area curve. ba_func <- bearing_area(normforest) # rotate the values and re-plot xval <- environment(ba_func)$y yval <- (1 - environment(ba_func)$x) plot(yval ~ xval)
Finds the best fit polynomial surface for a raster or matrix. This function tests least squares polynomial fits with orders of 0 - 3 and determines which order minimizes the error when the fit is subtracted from the original image.
bestfitplane(x)bestfitplane(x)
x |
A raster or matrix. |
A raster or matrix of the same size as the input with values predicted from the best polynomial fit.
# import raster image data(orforest) orforest <- terra::unwrap(orforest) # find the least squares polynomial surface poly <- bestfitplane(orforest) # plot the fit terra::plot(poly)# import raster image data(orforest) orforest <- terra::unwrap(orforest) # find the least squares polynomial surface poly <- bestfitplane(orforest) # plot the fit terra::plot(poly)
This function serves to shift the zero-frequency component of the Fourier transform to the center of the matrix.
fftshift(x, dim = -1)fftshift(x, dim = -1)
x |
An n x n Fourier transform matrix. |
dim |
Which dimension to shift the matrix. -1 swaps up/down and left/right. 1 swaps up/down. 2 swaps left/right. |
An n x n matrix with the zero-frequency component of the Fourier transform in the center.
#' This function was created from code posted by rayryeng at: https://stackoverflow.com/questions/38230794/how-to-write-fftshift-and-ifftshift-in-r.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # convert to matrix form M <- ncol(normforest) N <- nrow(normforest) zmat <- matrix(terra::values(normforest), ncol = M, nrow = N, byrow = TRUE) # calculate fourier transform and shift ftmat <- fft(zmat) ftshift <- fftshift(ftmat) # plot real component r <- terra::setValues(normforest, Re(ftshift)) terra::plot(r)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # convert to matrix form M <- ncol(normforest) N <- nrow(normforest) zmat <- matrix(terra::values(normforest), ncol = M, nrow = N, byrow = TRUE) # calculate fourier transform and shift ftmat <- fft(zmat) ftshift <- fftshift(ftmat) # plot real component r <- terra::setValues(normforest, Re(ftshift)) terra::plot(r)
Locates the flattest x percentage of the Bearing Area curve. Meant to locate the flattest 40 percent of the Bearing Area curve as used in several roughness parameter calculations.
find_flat(x, perc = 0.4)find_flat(x, perc = 0.4)
x |
A raster or matrix. |
perc |
Numeric between 0 and 1. The percentage of the curve over which to fit the line. |
A list containing the equation for the best fit line, the predicted values from that line, the high and low y-intercept values for the intersection points of the line with the Bearing Area curve, and the high and low x-intercept values for the intersection points of the line with the Bearing Area curve.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # locate the flattest 40% of the bearing area curve line_data <- find_flat(normforest, perc = 0.4) # extract the equation of the line bf_line <- line_data[[1]]# import raster image data(normforest) normforest <- terra::unwrap(normforest) # locate the flattest 40% of the bearing area curve line_data <- find_flat(normforest, perc = 0.4) # extract the equation of the line bf_line <- line_data[[1]]
Locates local peaks on a raster or matrix. A peak is defined as any pixel where all 8 surrounding pixels have lower values, and the center pixel has a positive value.
findpeaks(x)findpeaks(x)
x |
A raster or matrix. |
A dataframe of local peak locations (x, y) and
values (val). The raster or matrix location index (ind),
row (row), and column (col) are also listed.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # locate peaks peaks <- findpeaks(normforest) # calculate the summit density (# peaks/area) N <- ncol(normforest) M <- nrow(normforest) Sds <- nrow(peaks) / ((N - 1) * (M - 1))# import raster image data(normforest) normforest <- terra::unwrap(normforest) # locate peaks peaks <- findpeaks(normforest) # calculate the summit density (# peaks/area) N <- ncol(normforest) M <- nrow(normforest) Sds <- nrow(peaks) / ((N - 1) * (M - 1))
Locates local valleys on a raster or matrix. A valley is defined as any pixel where all 8 surrounding pixels have higher values, and the center pixel has a negative value.
findvalleys(x)findvalleys(x)
x |
A raster or matrix. |
A dataframe of local valley locations (x, y) and
values (val). The raster or matrix location index (ind),
row (row), and column (col) are also listed.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # locate peaks and valleys peaks <- findpeaks(normforest) valleys <- findvalleys(normforest) # find top 5 peaks, valleys top_peaks <- peaks[order(peaks$val, decreasing = TRUE)[1:5],] bottom_valleys <- valleys[order(valleys$val)[1:5],] # calculate the ten-point height S10z <- (sum(top_peaks$val) + sum(abs(bottom_valleys$val))) / 5# import raster image data(normforest) normforest <- terra::unwrap(normforest) # locate peaks and valleys peaks <- findpeaks(normforest) valleys <- findvalleys(normforest) # find top 5 peaks, valleys top_peaks <- peaks[order(peaks$val, decreasing = TRUE)[1:5],] bottom_valleys <- valleys[order(valleys$val)[1:5],] # calculate the ten-point height S10z <- (sum(top_peaks$val) + sum(abs(bottom_valleys$val))) / 5
Fits a polynomial surface of order n to a raster
or matrix.
fitplane(x, order)fitplane(x, order)
x |
A raster or matrix. |
order |
Numeric. Indicates the polynomial order to be fit. |
A matrix with values predicted from the polynomial fit.
# import raster image data(orforest) orforest <- terra::unwrap(orforest) # find the 2nd order least squares polynomial surface polyfit <- fitplane(orforest, order = 2) # create raster of polyfit x <- terra::setValues(orforest, polyfit) # plot the fit terra::plot(x)# import raster image data(orforest) orforest <- terra::unwrap(orforest) # find the 2nd order least squares polynomial surface polyfit <- fitplane(orforest, order = 2) # create raster of polyfit x <- terra::setValues(orforest, polyfit) # plot the fit terra::plot(x)
Calculates the surface area of a flat raster or matrix with the same x, y bounds as the study surface.
flatsa(x)flatsa(x)
x |
A raster or matrix. |
This function scales both x and y to between 0 and 1. This is done because most satellite data have units where the x, y units do not equal the z units and because the flat surface area is usually compared to the actual surface area. Surface area is calculated over the sample area (N-1, M-1).
A numeric value representing the scaled surface area of a flattened surface with the same x, y bounds.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate flattened surface area flatsa(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate flattened surface area flatsa(normforest)
Calculates the various texture metrics over windows centered
on individual pixels. This creates a continuous surface of the
texture metric.
This function is a modified version of the window_lsm function from the
landscapemetrics package (Hesselbarth et al. 2019).
focal_metrics(x, window, metrics, progress, ...)focal_metrics(x, window, metrics, progress, ...)
x |
A raster or matrix. Image over which to apply focal window calculations. |
window |
Matrix. The focal window used to create the image. |
metrics |
List. List of metrics to apply. Function names must be strings. |
progress |
Logical. Display progress through metrics list? |
... |
Additional arguments for the metric functions. All applicable arguments will be applied to the entire list of metrics. |
Metrics available from geodiv package:
'sa': average surface roughness
'sq': root mean square roughness
's10z': ten-point height
'sdq': root mean square slope of surface, 2-point method
'sdq6': root mean square slope of surface, 7-point method
'sdr': surface area ratio
'sbi': surface bearing index
'sci': core fluid retention index
'ssk': skewness
'sku': kurtosis
'sds': summit density
'sfd': 3d fractal dimension
'srw': dominant radial wavelength, radial wavelength index, mean half wavelength
'std': angle of dominating texture, texture direction index
'svi': valley fluid retention index
'stxr': texture aspect ratio
'ssc': mean summit curvature
'sv': maximum valley depth
'sph': maximum peak height
'sk': core roughness depth
'smean': mean peak height
'svk': reduced valley depth
'spk': reduced peak height
'scl': correlation length
'sdc': bearing area curve height interval
A raster of the metric calculated in windows over the raster or matrix. If the input was a matrix, the function will return a raster with an extent of [0, 1, 0, 1].
Hesselbarth, M.H.K., Sciaini, M., With, K.A., Wiegand, K., Nowosad, J. 2019. landscapemetrics: an open-source R tool to calculate landscape metrics. - Ecography 42:1648-1657(ver. 0).
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # crop raster to smaller area x <- terra::crop(normforest, terra::ext(normforest[1:100, 1:100, drop = FALSE])) # get a surface of root mean square roughness sa_img <- focal_metrics(x = x, window = matrix(1, 5, 5), metrics = list('sa'), progress = TRUE) # plot the result terra::plot(sa_img$sa)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # crop raster to smaller area x <- terra::crop(normforest, terra::ext(normforest[1:100, 1:100, drop = FALSE])) # get a surface of root mean square roughness sa_img <- focal_metrics(x = x, window = matrix(1, 5, 5), metrics = list('sa'), progress = TRUE) # plot the result terra::plot(sa_img$sa)
Determines the value of the bearing area curve for a
specific value along the x-axis (xval).
height_ba(x, xval)height_ba(x, xval)
x |
A raster or matrix. |
xval |
Numeric value along the x-axis. |
A numeric value of the bearing area function
corresponding to xval.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the bearing area function value # corresponding to an x value of 0.4 val <- height_ba(normforest, 0.4)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the bearing area function value # corresponding to an x value of 0.4 val <- height_ba(normforest, 0.4)
A raster image of Normalized Difference
Vegetation Index (NDVI) errors derived from Landsat
data for a small portion of SW Oregon state. This
raster was calculated by subtracting the best-fit polynomial
plane from the orforest values. The best-fit
polynomial plane was calculated using bestfitplane.
normforestnormforest
A raster image with 371 x 371 pixels
-0.5854638 – 0.1585918
-123, -122.9, 43.0002, 43.1
30m x 30m (0.002694946 degrees)
WGS84
1
...
NDVI values are derived from Landsat scene path 45, row 30 summarized as the mean NDVI value between June and August 2000 at roughly 30m resolution. Clouds were removed from the Landsat scene before calculating the mean. The image was created using Google Earth Engine in August 2018.
A raster image of Shuttle Radar Topography Mission (SRTM) elevation for a portion of southwestern Oregon.
orelevationorelevation
A raster image with 371 x 371 pixels
433 – 1390
-123.0001, -122.9002, 43.00015, 43.10013
30m x 30m (0.002694946 degrees)
WGS84
1
...
Elevation values are from the SRTM data for 2000 and are at roughly 30m resolution. The image was created using Google Earth Engine in October 2019.
A raster image of Normalized Difference Vegetation Index (NDVI) derived from Landsat data for a small portion of SW Oregon state.
orforestorforest
A raster image with 371 x 371 pixels
0 – 1
-123, -122.9, 43.0002, 43.1
30m x 30m (0.002694946 degrees)
WGS84
1
...
NDVI values are derived from Landsat scene path 45, row 30 summarized as the mean NDVI value between June and August 2000 at roughly 30m resolution. Clouds were removed from the Landsat scene before calculating the mean. The image was created using Google Earth Engine in August 2018.
Extends edge values of a raster or matrix by a specified number of pixels.
pad_edges(x, size, val = NULL)pad_edges(x, size, val = NULL)
x |
A matrix. |
size |
Numeric. Number of pixels to add to each side. |
val |
Numeric. If NULL (default), this extends the edge values out. If not null, this value will be used for the extended cells. |
A raster with edges padded size number of pixels on each edge.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # crop raster to much smaller area x <- pad_edges(as.matrix(normforest), 3, val = NA)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # crop raster to much smaller area x <- pad_edges(as.matrix(normforest), 3, val = NA)
Calculates and plots the Bearing Area curve for a raster
or matrix using the bearing_area() function (with correctly
rotated results).
plot_ba_curve(x, divisions = FALSE)plot_ba_curve(x, divisions = FALSE)
x |
A raster or matrix. |
divisions |
Logical, defaults to |
If divisions = TRUE, the lines representing the
best fit line to the flattest 40 percent of the curve will be
shown, as well as both the x and y interception points
of that line.
Plots the Bearing Area curve.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # plot the bearing area curve plot_ba_curve(normforest, divisions = TRUE)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # plot the bearing area curve plot_ba_curve(normforest, divisions = TRUE)
Finds the best fit polynomial surface for a raster or matrix and subtracts it from the actual values. The output image has positive values where the actual values are higher than the surface and negative values where the actual value are lower than the surface.
remove_plane(x)remove_plane(x)
x |
A raster or matrix. |
A raster or matrix of the same size as the input with values equal to the difference between the original and bestfit plane.
# import raster image data(orforest) orforest <- terra::unwrap(orforest) # remove the least squares polynomial surface new_rast <- remove_plane(orforest) # plot terra::plot(new_rast)# import raster image data(orforest) orforest <- terra::unwrap(orforest) # remove the least squares polynomial surface new_rast <- remove_plane(orforest) # plot terra::plot(new_rast)
Rotates a matrix 180 degrees. Code is from https://stackoverflow.com/questions/16496210/rotate-a-matrix-in-r-by-90-degrees-clockwise.
rotate(x)rotate(x)
x |
A matrix. |
A matrix rotate 180 degrees.
Calculates the average height above the mean surface for the five highest local maxima plus the average height below the mean surface for the five lowest local minima.
s10z(x)s10z(x)
x |
A raster or matrix. |
A numeric value representing the ten-point height.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate ten-point height. S10z <- s10z(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate ten-point height. S10z <- s10z(normforest)
Finds the average roughness of a surface (Sa) as the absolute deviation of surface heights from the mean surface height. Height is measured as the value of a raster and may not necessarily represent actual height.
sa(x)sa(x)
x |
A raster or matrix. |
A value of average roughness in the units of the original raster or matrix.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the surface roughness roughness <- sa(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the surface roughness roughness <- sa(normforest)
Determines the surface bearing index (Sbi), calculated as the ratio of root mean square roughness (Sq) to height at 5% of bearing area (z05).
sbi(x)sbi(x)
x |
A raster or matrix. |
A numeric value representing the surface bearing index.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the surface bearing index Sbi <- sbi(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the surface bearing index Sbi <- sbi(normforest)
Determines the core fluid retention index (Sci). This value is the void volume (area under the bearing area curve) in the 'core' zone. See Figure 2a from Kedron et al. (2018) for more details.
sci(x)sci(x)
x |
A raster or matrix. |
A numeric value representing the core fluid retention index.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the core fluid retention index Sci <- sci(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the core fluid retention index Sci <- sci(normforest)
Calculates the smallest and largest distances to specified autocorrelation values (e.g., 0.2) of the areal autocorrelation function (AACF). All 180 degrees from the origin of the AACF image are considered for the calculation.
scl(x, threshold = c(0.2, 1/exp(1)), create_plot = FALSE)scl(x, threshold = c(0.2, 1/exp(1)), create_plot = FALSE)
x |
A raster or matrix. |
threshold |
A numeric vector containing values between 0 and 1. Indicates the autocorrelation values to which the rates of decline are measured. |
create_plot |
Logical. Defaults to |
A list containing the minimum and maximum distances from an
autocorrelation value of 1 to the specified autocorrelation values < 1.
Distances are in the units of the x, y coordinates of the raster image. If more
than one threshold value is specified, the order of this list will be
[minval(t1), minval(t2), maxval(t1), maxval(t2)].
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate Scl20, the minimum distance to an autocorrelation value of 0.2 in the AACF Scl20 <- scl(normforest)[1]# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate Scl20, the minimum distance to an autocorrelation value of 0.2 in the AACF Scl20 <- scl(normforest)[1]
Determines the height interval (height distance) for points along the bearing area curve as defined by their x values.
sdc(x, low, high)sdc(x, low, high)
x |
A raster or matrix. |
low |
Numeric value along the x-axis corresponding to the lowest value of interest along the x-axis. |
high |
Numeric value along the y-axis corresponding to the highest value of interest along the x-axis. |
A numeric value of the difference in height of the y values along the bearing area curve corresponding to the specified x values.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the 10-40% height interval of the # bearing area curve val <- sdc(normforest, 0.1, 0.4)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the 10-40% height interval of the # bearing area curve val <- sdc(normforest, 0.1, 0.4)
Calculates the root mean square slope of a raster or matrix surface using the two-point method.
sdq(x)sdq(x)
x |
A raster or matrix. |
A numeric value representing the two-point root mean square slope, Sdq. The units of the returned value are change in z per one unit (pixel).
This function is based on the equations found at https://www.ntmdt-si.ru/data/media/files/manuals/image_analisys_p9_nov12.e.pdf.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate root mean square slope Sdq <- sdq(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate root mean square slope Sdq <- sdq(normforest)
Calculates the area root mean square slope of a raster or matrix surface using the seven-point method.
sdq6(x)sdq6(x)
x |
A raster or matrix. |
A numeric value representing the seven-point root mean square slope, Sdq6. The units of the returned value are change in z per one unit (pixel).
This function is based on the equations found at https://www.ntmdt-si.ru/data/media/files/manuals/image_analisys_p9_nov12.e.pdf.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate area root mean square slope Sdq6 <- sdq6(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate area root mean square slope Sdq6 <- sdq6(normforest)
Calculates the surface area ratio of a raster or matrix. This is the ratio of a flat surface to the actual surface.
sdr(x)sdr(x)
x |
A raster or matrix. |
This function scales both x and y, as well as the surface value (z), to between 0 and 1 to best match their units. This is done because most satellite data have units where the x, y units do not equal the z units. Surface area is calculated over the sample area (N-1, M-1).
A numeric value representing the surface area ratio.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate the surface area ratio Sdr <- sdr(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate the surface area ratio Sdr <- sdr(normforest)
Calculates the summit density of a raster or matrix. Summit density is the number of local peaks per unit area.
sds(x)sds(x)
x |
A raster or matrix. |
A numeric value representing the summit density.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate summit density. Sds <- sds(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate summit density. Sds <- sds(normforest)
Calculates the 3D fractal dimension of a raster using the triangular prism surface area method.
sfd(x, silent = FALSE)sfd(x, silent = FALSE)
x |
A raster or matrix. |
silent |
Logical. If |
A numeric value representing the fractal dimension of the image.
Clarke, K.C., 1986. Computation of the fractal dimension of topographic surfaces using the triangular prism surface area method. Computers & Geosciences, 12(5), pp.713-722.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate the fractal dimension Sfd <- sfd(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate the fractal dimension Sfd <- sfd(normforest)
Calculates the 3D fractal dimension of a raster using the triangular prism surface area method.
sfd_(mat)sfd_(mat)
mat |
A matrix. |
A numeric value representing the fractal dimension of the image.
Clarke, K.C., 1986. Computation of the fractal dimension of topographic surfaces using the triangular prism surface area method. Computers & Geosciences, 12(5), pp.713-722.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # convert to matrix mat <- matrix(normforest[], ncol = ncol(normforest), nrow = nrow(normforest)) # calculate the fractal dimension Sfd <- sfd_(mat)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # convert to matrix mat <- matrix(normforest[], ncol = ncol(normforest), nrow = nrow(normforest)) # calculate the fractal dimension Sfd <- sfd_(mat)
Calculates the area below a curve from
points a to b. This function is provided
for general use.
simpsons(f, a, b, n = 100)simpsons(f, a, b, n = 100)
f |
A function. |
a |
Numeric. The left x boundary. |
b |
Numeric. The right x boundary. |
n |
Numeric. The number of subdivisions along the function line. |
Note that if y-values are negative, this returns the area above the function line.
A numeric value representing the area under the curve with
x bounds a and b.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # basic values z <- terra::values(normforest) # calculate cumulative probability density function of surface 'height' (= ndvi) mod <- ecdf((1 - z)) # calculate integral int_area <- simpsons(f = mod, b = 1, a = 0.8, n = 500)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # basic values z <- terra::values(normforest) # calculate cumulative probability density function of surface 'height' (= ndvi) mod <- ecdf((1 - z)) # calculate integral int_area <- simpsons(f = mod, b = 1, a = 0.8, n = 500)
Determines the core roughness depth (Sk), the height difference between y values of the intersection points of the least mean square line fit to the flattest 40% of the bearing area curve. See Figure 2a from Kedron et al. (2018) for more details.
sk(x)sk(x)
x |
A raster. |
A numeric value representing the core roughness depth of the image.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the core roughness depth Sk <- sk(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the core roughness depth Sk <- sk(normforest)
Finds the kurtosis for a distribution of raster or matrix values (Sku). Kurtosis represents the peakedness of the raster surface height distribution. Height is measured as the value of a raster/matrix and may not necessarily represent actual height.
sku(x, excess = TRUE)sku(x, excess = TRUE)
x |
A raster or matrix. |
excess |
Logical, defaults to |
A numeric value representing kurtosis.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the excess kurtosis of the raster distribution Sku <- sku(normforest, excess = TRUE)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the excess kurtosis of the raster distribution Sku <- sku(normforest, excess = TRUE)
Calculates the slopes along the bearing area curve of a raster or matrix. Slopes are determined at points x, from point x - h to x + h.
slopecalc(x, h, f)slopecalc(x, h, f)
x |
A vector of x values. |
h |
Spacing before and after each point. 2h is the distance over which slopes are calculated. |
f |
Bearing area function as calculated with bearing_area. |
A dataframe with the slope for each segment with centerpoint x.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the slopes along the bearing area curve ba <- bearing_area(normforest) x <- seq(0, 1, length.out = 100000) slopes <- slopecalc(x = x, h = 0.01, f = ba)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the slopes along the bearing area curve ba <- bearing_area(normforest) x <- seq(0, 1, length.out = 100000) slopes <- slopecalc(x = x, h = 0.01, f = ba)
Calculates the average slope over every segment of a specified percentage length of the total bearing area curve.
slopemeans(slopes, l = 0.4)slopemeans(slopes, l = 0.4)
slopes |
A dataframe containing all slopes along the bearing area curve, calculated using the slopecalc function. |
l |
Percentage of the curve over which to calculate mean slope. |
A dataframe with the average slope over segments beginning at specified x locations along the bearing area curve. 'slope' represents the mean slope over the segment, 'xstart' is the beginning x location of the segment, and 'xend' is the concluding x location of the segment.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the average slope of segments of the bearing area # curve. ba <- bearing_area(normforest) x <- seq(0, 1, length.out = 10000) slopes <- slopecalc(x = x, h = 0.01, f = ba) slopes_forty <- slopemeans(slopes = slopes, l = 0.4)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the average slope of segments of the bearing area # curve. ba <- bearing_area(normforest) x <- seq(0, 1, length.out = 10000) slopes <- slopecalc(x = x, h = 0.01, f = ba) slopes_forty <- slopemeans(slopes = slopes, l = 0.4)
Finds the mean height of positive values in the landscape (mean peak height; Smean) for a raster or matrix representing a surface.
smean(x)smean(x)
x |
A raster or matrix. |
A numeric value of mean peak height.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the maximum peak height Smean <- smean(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the maximum peak height Smean <- smean(normforest)
Finds the absolute value of the highest value in the landscape (maximum peak height; Sph) for a raster or matrix representing a surface.
sph(x)sph(x)
x |
A raster or matrix. |
A numeric value of maximum peak height.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the maximum peak height Sph <- sph(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the maximum peak height Sph <- sph(normforest)
Determines the reduced peak height (Spk), the height difference between the maximum y value of the bearing area curve and the y value of the highest intersection point of the least mean square line fit to the flattest 40% of the bearing area curve. See Figure 2a from Kedron et al. (2018) for more details.
spk(x)spk(x)
x |
A raster or matrix. |
A numeric value representing the reduced peak height.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the reduced peak height Spk <- spk(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the reduced peak height Spk <- spk(normforest)
Finds the root mean square roughness of a surface (Sq) as the standard deviation of surface heights from the mean surface height. Height is measured as the value of a raster and may not necessarily represent actual height.
sq(x)sq(x)
x |
A raster or matrix. |
A value of root mean square roughness in the units of the original raster or matrix.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the surface roughness roughness <- sq(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the surface roughness roughness <- sq(normforest)
Calculates the dominant radial wavelength, radial wavelength index, and mean half wavelength of the radial Fourier spectrum. See Kedron et al. (2018) for more detailed description.
srw(x, create_plot = FALSE, option = c(1, 2, 3))srw(x, create_plot = FALSE, option = c(1, 2, 3))
x |
A raster or matrix. |
create_plot |
Logical. If |
option |
Numeric. Code for which output metric(s) to return. 1 = Srw, 2 = Srwi, 3 = Shw. |
A vector containing numeric values for the dominant radial wavelength, radial wavelength index, and mean half wavelength.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate metrics srwvals <- srw(normforest) # extract each value Srw <- srwvals[1] Srwi <- srwvals[2] Shw <- srwvals[3]# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate metrics srwvals <- srw(normforest) # extract each value Srw <- srwvals[1] Srwi <- srwvals[2] Shw <- srwvals[3]
Calculates the mean summit curvature of a raster or matrix. Mean summit curvature is the average principle curvature of local maximas on the surface.
ssc(x)ssc(x)
x |
A raster or matrix. |
A numeric value representing the average curvature of surface peaks.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate mean summit curvature Ssc <- ssc(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate mean summit curvature Ssc <- ssc(normforest)
Finds the Fisher-Pearson coefficient of skewness for raster or matrix values (Ssk). Skewness represents the asymmetry of the surface height distribution. Height is measured as the value of a raster/matrix and may not necessarily represent actual height.
ssk(x, adj = TRUE)ssk(x, adj = TRUE)
x |
A raster or matrix. |
adj |
Logical, defaults to |
A numeric value representing skewness. # import raster image data(normforest) normforest <- terra::unwrap(normforest)
# find the adjusted coefficient of skewness Ssk <- ssk(normforest, adj = TRUE)
Calculates the angle of dominating texture and the texture direction index of the Fourier spectrum image calculated from a raster image (see Kedron et al. 2018).
std(x, create_plot = FALSE, option = c(1, 2))std(x, create_plot = FALSE, option = c(1, 2))
x |
A raster or matrix. |
create_plot |
Logical. If |
option |
Numeric. Code for which output metric(s) to return. 1 = Std, 2 = Stdi. |
A vector containing numeric values for the angle of dominating texture and the texture direction index.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate Std and Stdi stdvals <- std(normforest) # extract each value Std <- stdvals[1] Stdi <- stdvals[2]# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate Std and Stdi stdvals <- std(normforest) # extract each value Std <- stdvals[1] Stdi <- stdvals[2]
Calculates the texture aspect ratio (Str) at defined autocorrelation values. The texture aspect ratio is the ratio of the fastest to the slowest decay lengths of the autocorrelation function to the defined autocorrelation values.
stxr(x, threshold = c(0.2, 1/exp(1)))stxr(x, threshold = c(0.2, 1/exp(1)))
x |
A raster or matrix. |
threshold |
A vector of autocorrelation values with values between 0 and 1. Indicates the autocorrelation value(s) to which the rates of decline are measured. |
A vector with length equal to that of threshold
containing the texture aspect ratio(s) for the input autocorrelation
value(s).
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # estimate the texture aspect ratio for autocorrelation # thresholds of 0.20 and 0.37 (1/e) strvals <- stxr(normforest, threshold = c(0.20, 1 / exp(1))) # calculate Str20, the texture aspect ratio for # autocorrelation value of 0.2 in the AACF Str20 <- strvals[1]# import raster image data(normforest) normforest <- terra::unwrap(normforest) # estimate the texture aspect ratio for autocorrelation # thresholds of 0.20 and 0.37 (1/e) strvals <- stxr(normforest, threshold = c(0.20, 1 / exp(1))) # calculate Str20, the texture aspect ratio for # autocorrelation value of 0.2 in the AACF Str20 <- strvals[1]
Calculates the scaled surface area of a raster or matrix.
surface_area(x)surface_area(x)
x |
A raster or matrix. |
This function scales both x and y, as well as the surface value (z), to between 0 and 1 to best match their units. This is done because most satellite data have units where the x, y units do not equal the z units. The surface area represents the surface area of the sample area (N-1, M-1).
Note that the surface object may have NA values around the edges, but should not have any missing values within the main area.
A numeric value representing the scaled surface area.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate surface area surface_area(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # calculate surface area surface_area(normforest)
Finds the absolute value of the lowest value in the landscape (maximum valley depth; Sv) for a raster or matrix representing a surface.
sv(x)sv(x)
x |
A raster or matrix. |
A numeric value of maximum valley depth.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the maximum valley depth Sv <- sv(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # find the maximum valley depth Sv <- sv(normforest)
Determines the valley fluid retention index (Svi). This value is the void volume (area under the bearing area curve) in the 'valley' zone. See Figure 2a from Kedron et al. (2018) for more details.
svi(x)svi(x)
x |
A raster or matrix. |
A numeric value representing the valley fluid retention index.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the valley fluid retention index Svi <- svi(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the valley fluid retention index Svi <- svi(normforest)
Determines the reduced valley depth (Svk), the height difference between y value of the lowest intersection point of the least mean square line fit to the flattest 40% of the bearing area curve and the minimum y value of the bearing area curve. See Figure 2a from Kedron et al. (2018) for more details.
svk(x)svk(x)
x |
A raster or matrix. |
A numeric value representing the reduced valley depth.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the reduced valley depth Svk <- svk(normforest)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # determine the reduced valley depth Svk <- svk(normforest)
Calculates the various texture metrics over windows centered on individual pixels. This creates a continuous surface of the texture metric.
texture_image( x, window_type = "square", size = 5, in_meters = FALSE, metric, args = NULL, parallel = TRUE, ncores = NULL, nclumps = 100 )texture_image( x, window_type = "square", size = 5, in_meters = FALSE, metric, args = NULL, parallel = TRUE, ncores = NULL, nclumps = 100 )
x |
A raster or matrix. |
window_type |
Character. Type of window, either circular or square. |
size |
Numeric. Size of window (edge length) or diameter (in meters). |
in_meters |
Logical. Is the size given in meters? |
metric |
Character. Metric to calculate for each window. Metrics from the geodiv package are listed below. |
args |
List. Arguments from function to be applied over each window (e.g., list(threshold = 0.2)). |
parallel |
Logical. Option to run the calculations in parallel on available cores. |
ncores |
Numeric. If parallel is TRUE, number of cores on which to run the calculations. Defaults to all available, minus 1. |
nclumps |
Numeric. Number of clumps to split the raster or matrix into. |
Note that this function is meant to work on rasters with an equal area projection.
Metrics available from geodiv package:
'sa': average surface roughness
'sq': root mean square roughness
's10z': ten-point height
'sdq': root mean square slope of surface, 2-point method
'sdq6': root mean square slope of surface, 7-point method
'sdr': surface area ratio
'sbi': surface bearing index
'sci': core fluid retention index
'ssk': skewness
'sku': kurtosis
'sds': summit density
'sfd': 3d fractal dimension
'srw': dominant radial wavelength, radial wavelength index, mean half wavelength
'std': angle of dominating texture, texture direction index
'svi': valley fluid retention index
'stxr': texture aspect ratio
'ssc': mean summit curvature
'sv': maximum valley depth
'sph': maximum peak height
'sk': core roughness depth
'smean': mean peak height
'svk': reduced valley depth
'spk': reduced peak height
'scl': correlation length
'sdc': bearing area curve height interval
A raster or list of rasters (if function results in multiple outputs) with pixel values representative of the metric value for the window surrounding that pixel.
The total window size for square windows will be (size * 2) + 1.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # crop raster to smaller area x <- terra::crop(normforest, terra::ext(normforest[1:100, 1:100, drop = FALSE])) # get a surface of root mean square roughness sa_img <- texture_image(x = x, window = 'square', size = 5, metric = 'sa', parallel = TRUE, ncores = 1, nclumps = 20) # plot the result terra::plot(sa_img)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # crop raster to smaller area x <- terra::crop(normforest, terra::ext(normforest[1:100, 1:100, drop = FALSE])) # get a surface of root mean square roughness sa_img <- texture_image(x = x, window = 'square', size = 5, metric = 'sa', parallel = TRUE, ncores = 1, nclumps = 20) # plot the result terra::plot(sa_img)
Calculates the various texture metrics over a window centered on an individual pixel.
window_metric( x, coords, window_type = "square", size = 11, metric, args = NULL )window_metric( x, coords, window_type = "square", size = 11, metric, args = NULL )
x |
A raster or matrix. |
coords |
Dataframe. Coordinates of window edges. |
window_type |
Character. Type of window, either circular or square. |
size |
Numeric. Edge length or diameter of window in number of pixels. |
metric |
Character. Metric to calculate for each window. Metrics from the geodiv package are listed below. |
args |
List. Arguments from function to be applied over each window (e.g., list(threshold = 0.2)). |
Metrics from geodiv package:
'sa': average surface roughness
'sq': root mean square roughness
's10z': ten-point height
'sdq': root mean square slope of surface, 2-point method
'sdq6': root mean square slope of surface, 7-point method
'sdr': surface area ratio
'sbi': surface bearing index
'sci': core fluid retention index
'ssk': skewness
'sku': kurtosis
'sds': summit density
'sfd': 3d fractal dimension
'srw': dominant radial wavelength, radial wavelength index, mean half wavelength
'std': angle of dominating texture, texture direction index
'svi': valley fluid retention index
'stxr': texture aspect ratio
'ssc': mean summit curvature
'sv': maximum valley depth
'sph': maximum peak height
'sk': core roughness depth
'smean': mean peak height
'svk': reduced valley depth
'spk': reduced peak height
'scl': correlation length
'sdc': bearing area curve height interval
A raster with pixel values representative of the metric value for the window surrounding that pixel.
Note that if calculating the metric at the edge of a raster or matrix,
the input raster/matrix must be padded. This can be done using the pad_edges
function.
Calculates a matrix of values with a negative or positive, x or y, offset.
zshift(r, xdist = 0, ydist = 0, xrm, yrm, scale = FALSE)zshift(r, xdist = 0, ydist = 0, xrm, yrm, scale = FALSE)
r |
A raster or matrix. |
xdist |
Numeric indicating the number and direction (+, -) of columns for the offset. |
ydist |
Numeric indicating the number and direction (+, -) of rows for the offset. |
xrm |
Numeric value or vector indicating the number of
columns to be removed from the final matrix. If not set,
this value defaults to |
yrm |
Numeric value or vector indicating the number
of rows to be removed from the final matrix. If not set,
this value defaults to |
scale |
Logical. Indicates whether or not to scale the values of the raster. |
A numeric vector of values created from a matrix of the values
with the specified offset. The vector is created from a matrix with
xrm fewer columns and yrm fewer rows than the original
raster value matrix.
# import raster image data(normforest) normforest <- terra::unwrap(normforest) # remove right and bottom borders 2 deep noborder <- zshift(normforest, xdist = 2, ydist = 2)# import raster image data(normforest) normforest <- terra::unwrap(normforest) # remove right and bottom borders 2 deep noborder <- zshift(normforest, xdist = 2, ydist = 2)