Skip to contents

Calculates surface area (Jenness, 2004) to planar area rugosity and by default corrects planar area for slope using the arc-chord ratio (Du Preez, 2015). Additionally, the method has been modified to allow for calculations at multiple different window sizes (see details and Ilich et al. (2023)).

Usage

SAPA(
  r = NULL,
  w = 1,
  slope_correction = TRUE,
  na.rm = FALSE,
  include_scale = FALSE,
  slope_layer = NULL,
  sa_layer = NULL,
  filename = NULL,
  overwrite = FALSE,
  wopt = list()
)

Arguments

r

DTM as a SpatRaster or RasterLayer in a projected coordinate system where map units match elevation/depth units (Not used if both slope_layer and sa_layer are specified).

w

A single number or a vector of length 2 (row, column) specifying the dimensions of the rectangular window over which surface area will be summed. Window size must be an odd number. 1 refers to "native" scale and surface area and planar area will be calculated per cell (the traditional implementation).

slope_correction

Whether to use the arc-chord ratio to correct planar area for slope (default is TRUE)

na.rm

logical indicating whether to remove/account for NAs in calculations. If FALSE any calculations involving NA will be NA. If TRUE, NA values will be removed and accounted for.

include_scale

logical indicating whether to append window size to the layer names (default = FALSE)

slope_layer

Optionally specify an appropriate slope layer IN RADIANS to use. If not supplied, it will be calculated using the SlpAsp function based on Misiuk et al (2021). The slope layer should have a window size that is 2 larger than the w specified for SAPA.

sa_layer

Optionally specify a surface area raster that contains the surface area on a per cell level. This can be calculated with the SurfaceArea function. If calculating SAPA at multiple scales it will be more efficient to supply this so that it does not need to be calculated every time.

filename

character Output filename.

overwrite

logical. If TRUE, filename is overwritten (default is FALSE).

wopt

list with named options for writing files as in writeRaster

Value

a SpatRaster or RasterLayer

Details

Planar area is calculated as the x_dis * y_dis if uncorrected for slope and (x_dis * y_dis)/cos(slope) if corrected for slope. When w=1, this is called "native" scale and is equivalent to what is presented in Du Preez (2015) and available in the ArcGIS Benthic Terrain Modeller add-on. In this case operations are performed on a per cell basis where x_dis is the resolution of the raster in the x direction (left/right) and y_dis is the resolution of the raster in the y direction (up/down) and slope is calculated using the Horn (1981) method. To expand this to multiple scales of analysis, at w > 1 slope is calculated based on Misiuk et al (2021) which provides a modification of the Horn method to extend the matric to multiple spatial scales. Planar area is calculated the same way as for w=1 except that now x_dis is the x resolution of the raster * the number of columns in the focal window, and y_dis is y resolution of the raster * the number of rows. For w > 1, surface area is calculated as the sum of surface areas within the focal window. Although the (modified) Horn slope is used by default to be consistent with Du Preez (2015), slope calculated using a different algorithm (e.g. Wood 1996) could be supplied using the slope_layer argument. Additionally, a slope raster can be supplied if you have already calculated it and do not wish to recalculate it. However, be careful to supply a slope layer measured in radians and calculated at the relevant scale (2 larger than the w of SAPA).

References

Du Preez, C., 2015. A new arc–chord ratio (ACR) rugosity index for quantifying three-dimensional landscape structural complexity. Landscape Ecol 30, 181–192. https://doi.org/10.1007/s10980-014-0118-8

Horn, B.K., 1981. Hill Shading and the Reflectance Map. Proceedings of the IEEE 69, 14-47.

Ilich, A. R., Misiuk, B., Lecours, V., & Murawski, S. A. (2023). MultiscaleDTM: An open-source R package for multiscale geomorphometric analysis. Transactions in GIS, 27(4). https://doi.org/10.1111/tgis.13067

Jenness, J.S., 2004. Calculating landscape surface area from digital elevation models. Wildlife Society Bulletin 32, 829-839.

Misiuk, B., Lecours, V., Dolan, M.F.J., Robert, K., 2021. Evaluating the Suitability of Multi-Scale Terrain Attribute Calculation Approaches for Seabed Mapping Applications. Marine Geodesy 44, 327-385. https://doi.org/10.1080/01490419.2021.1925789

Examples

r<- rast(volcano, extent= ext(2667400, 2667400 + 
ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), 
crs = "EPSG:27200")
sapa<- SAPA(r, w=c(5,5), slope_correction = TRUE)
plot(sapa)