SAPA Updates

- Add na.rm to SurfaceArea and SAPA by scaling mean to number of triangles/cells
- Also in C_AdjSD_narmT remove na_rm argument since it is not used

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: MultiscaleDTM
 Title: Multi-Scale Geomorphometric Terrain Attributes
-Version: 0.8.2
+Version: 0.8.3
 Authors@R: 
      c(
     person(given = "Alexander",

R/AdjSD.R

@@ -70,7 +70,7 @@ AdjSD<- function(r, w=c(3,3), na.rm=FALSE, include_scale=FALSE, filename=NULL, o
 
   #Fit Quadratic and Extract Residuals
   if(na.rm){
-    out<- terra::focalCpp(r, w=w, fun = C_AdjSD_narmT, X_full= X, na_rm=TRUE, fillvalue=NA, wopt=wopt)
+    out<- terra::focalCpp(r, w=w, fun = C_AdjSD_narmT, X_full= X, fillvalue=NA, wopt=wopt)
   } else{
     out<- terra::focalCpp(r, w=w, fun = C_AdjSD_narmF, X= X, Xt=Xt, XtX_inv=XtX_inv, fillvalue=NA, wopt=wopt)
     }

R/RcppExports.R

@@ -37,8 +37,8 @@ C_Qfit2_narmF <- function(z, X, Xt, XtX_inv, ni, nw) {
     .Call(`_MultiscaleDTM_C_Qfit2_narmF`, z, X, Xt, XtX_inv, ni, nw)
 }
 
-C_AdjSD_narmT <- function(z, X_full, na_rm, ni, nw) {
-    .Call(`_MultiscaleDTM_C_AdjSD_narmT`, z, X_full, na_rm, ni, nw)
+C_AdjSD_narmT <- function(z, X_full, ni, nw) {
+    .Call(`_MultiscaleDTM_C_AdjSD_narmT`, z, X_full, ni, nw)
 }
 
 C_AdjSD_narmF <- function(z, X, Xt, XtX_inv, ni, nw) {
@@ -49,8 +49,8 @@ C_TriArea <- function(a, b, c) {
     .Call(`_MultiscaleDTM_C_TriArea`, a, b, c)
 }
 
-C_SurfaceArea <- function(z, x_res, y_res, ni, nw) {
-    .Call(`_MultiscaleDTM_C_SurfaceArea`, z, x_res, y_res, ni, nw)
+C_SurfaceArea <- function(z, x_res, y_res, na_rm, ni, nw) {
+    .Call(`_MultiscaleDTM_C_SurfaceArea`, z, x_res, y_res, na_rm, ni, nw)
 }
 
 C_CountVals <- function(z, ni, nw) {

R/SAPA.R

@@ -1,11 +1,13 @@
 #' Calculates surface area to planar area rugosity
 #'
 #' 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)).
-#' @param r DTM as a SpatRaster or RasterLayer in a projected coordinate system where map units match elevation/depth units
+#' @param 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).
 #' @param 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).
 #' @param slope_correction Whether to use the arc-chord ratio to correct planar area for slope (default is TRUE)
+#' @param 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.
 #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
 #' @param 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.
+#' @param 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.
 #' @param filename character Output filename.
 #' @param overwrite logical. If TRUE, filename is overwritten (default is FALSE).
 #' @param wopt list with named options for writing files as in writeRaster
@@ -15,7 +17,7 @@
 #' crs = "EPSG:27200")
 #' sapa<- SAPA(r, w=c(5,5), slope_correction = TRUE)
 #' plot(sapa)
-#' @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). 
+#' @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).
 #' @return a SpatRaster or RasterLayer
 #' @import terra
 #' @importFrom raster raster
@@ -32,24 +34,68 @@
 #' 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
 #' @export
 
-SAPA<- function(r, w = 1, slope_correction=TRUE, include_scale=FALSE, slope_layer= NULL, filename=NULL, overwrite=FALSE, wopt=list()){
-  og_class<- class(r)[1]
-  if(og_class=="RasterLayer"){
-    r<- terra::rast(r) #Convert to SpatRaster
+SAPA<- function(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()){
+  if(is.null(r) & (is.null(slope_layer) | is.null(sa_layer))){
+    stop("Error: If r is NULL then 'slope_layer' and 'sa_layer' must be specified")
   }
+
+  og_class<- c(class(r)[1],class(slope_layer)[1], class(sa_layer)[1])
+
+  #Convert to SpatRaster
+  if(og_class[1]=="RasterLayer" & og_class[1]!="NULL"){
+    r<- terra::rast(r)
+  }
+  if(og_class[2]=="RasterLayer" & og_class[2]!="NULL"){
+    slope_layer<- terra::rast(slope_layer)
+  }
+  if(og_class[3]=="RasterLayer" & og_class[3]!="NULL"){
+    sa_layer<- terra::rast(sa_layer)
+  }
+
   #Input checks
-  if(!(og_class %in% c("RasterLayer", "SpatRaster"))){
+  if(any(!(og_class %in% c("RasterLayer", "SpatRaster", "NULL")))){
     stop("Error: Input must be a 'SpatRaster' or 'RasterLayer'")
   }
-  if(terra::nlyr(r)!=1){
-    stop("Error: Input raster must be one layer.")
+
+  if(!is.null(r)){
+    if(terra::nlyr(r)!=1){
+      stop("Error: Input raster 'r' must be one layer.")
+    }
+    if(isTRUE(terra::is.lonlat(r, perhaps=FALSE))){
+      stop("Error: Coordinate system of 'r' is Lat/Lon. Coordinate system must be projected with elevation/depth units matching map units.")
+    }
+    if(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){
+      warning("Coordinate system of 'r' may be Lat/Lon. Please ensure that the coordinate system is projected with elevation/depth units matching map units.")
+    }
   }
-  if(isTRUE(terra::is.lonlat(r, perhaps=FALSE))){
-    stop("Error: Coordinate system is Lat/Lon. Coordinate system must be projected with elevation/depth units matching map units.")
+
+  if(!is.null(slope_layer)){
+    if(terra::nlyr(slope_layer)!=1){
+      stop("Error: Input raster 'slope_layer' must be one layer.")
+    }
+    if(isTRUE(terra::is.lonlat(slope_layer, perhaps=FALSE))){
+      stop("Error: Coordinate system of 'slope_layer' is Lat/Lon. Coordinate system must be projected.")
+    }
+    if(terra::is.lonlat(slope_layer, perhaps=TRUE, warn=FALSE)){
+      warning("Coordinate system of 'slope_layer' may be Lat/Lon. Please ensure that the coordinate system is projected.")
+    }
+    if(terra::global(slope_layer, max, na.rm=TRUE) > 1.6){
+      stop("Error: Slope must be in radians")
+    }
   }
-  if(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){
-    warning("Coordinate system may be Lat/Lon. Please ensure that the coordinate system is projected with elevation/depth units matching map units.")
+
+  if(!is.null(sa_layer)){
+    if(terra::nlyr(sa_layer)!=1){
+      stop("Error: Input raster 'sa_layer' must be one layer.")
+    }
+    if(isTRUE(terra::is.lonlat(sa_layer, perhaps=FALSE))){
+      stop("Error: Coordinate system of 'sa_layer' is Lat/Lon. Coordinate system must be projected with units corresponding to map units^2.")
+    }
+    if(terra::is.lonlat(sa_layer, perhaps=TRUE, warn=FALSE)){
+      warning("Coordinate system of 'sa_layer' may be Lat/Lon. Please ensure that the coordinate system is projected with units corresponding to map units^2.")
+    }
   }
+
   if(length(w)==1){w<- rep(w,2)}
   if(length(w) > 2){
     stop("Specified window exceeds 2 dimensions")}
@@ -59,27 +105,34 @@ SAPA<- function(r, w = 1, slope_correction=TRUE, include_scale=FALSE, slope_laye
   if(all(w==c(1,1))){
     is_native<- TRUE} else{
       is_native<- FALSE} #Indicate whether SAPA is calculated at native scale
-
-  x_res<- terra::xres(r)
-  y_res<- terra::yres(r)
 
-  SA<- SurfaceArea(r, wopt=wopt)
+  if(is.null(sa_layer)){
+    sa_layer<- SurfaceArea(r, na.rm=na.rm, wopt=wopt)
+  }
 
   if(!is_native) {
-    SA<- terra::focal(SA, w= w, fun=sum, na.rm=FALSE, wopt=wopt)
-  } # Sum up surface area in focal window
+    if(na.rm){
+      sa_layer<- prod(w)*terra::focal(sa_layer, w= w, fun=mean, na.rm=TRUE, wopt=wopt)
+      } else{
+        sa_layer<- terra::focal(sa_layer, w= w, fun=sum, na.rm=FALSE, wopt=wopt)
+      }
+  } # Sum up surface area in focal window (or scale mean to number of cells)
 
   if(is.null(slope_layer) & slope_correction){
     if (all(w==1)){
       slope_layer<- terra::terrain(r, v="slope", unit="radians", neighbors=8, wopt=wopt)
       } else{
-        slope_layer<- SlpAsp(r, w=w+2, unit="radians", method="queen", metrics= "slope", na.rm=FALSE, include_scale=FALSE, wopt=wopt)
+        slope_layer<- SlpAsp(r, w=w+2, unit="radians", method="queen", metrics= "slope", na.rm=na.rm, include_scale=FALSE, wopt=wopt)
       }
-    }
-  PA<- (x_res*w[2]) * (y_res*w[1])
-  if(slope_correction){PA<- PA/terra::math(slope_layer, fun="cos", wopt=wopt)}#Planar area corrected for slope
+  }
+
+  x_res<- terra::xres(slope_layer)
+  y_res<- terra::yres(slope_layer)
+
+  pa_layer<- (x_res*w[2]) * (y_res*w[1])
+  if(slope_correction){pa_layer<- pa_layer/terra::math(slope_layer, fun="cos", wopt=wopt)}#Planar area corrected for slope
 
-  sapa<- SA/PA
+  sapa<- sa_layer/pa_layer
   names(sapa)<- "sapa"
 
   if(include_scale){
@@ -90,7 +143,7 @@ SAPA<- function(r, w = 1, slope_correction=TRUE, include_scale=FALSE, slope_laye
     }}
 
   #Return
-  if(og_class =="RasterLayer"){
+  if(!("SpatRaster" %in% og_class)){
     sapa<- raster::raster(sapa)
     if(!is.null(filename)){
       return(raster::writeRaster(sapa, filename=filename, overwrite=overwrite))

R/SurfaceArea.R

@@ -2,6 +2,7 @@
 #'
 #' Calculates surface area on a per cell basis of a DTM based on Jenness, 2004.
 #' @param r DTM as a SpatRaster or RasterLayer in a projected coordinate system where map units match elevation/depth units
+#' @param na.rm Logical indicating whether to remove NAs from calculations. When FALSE, the sum of the eight triangles is calculated. When TRUE, the mean of the created triangles is calculated and multiplied by 8 to scale it to the proper area.
 #' @param filename character Output filename.
 #' @param overwrite logical. If TRUE, filename is overwritten (default is FALSE).
 #' @param wopt list with named options for writing files as in writeRaster
@@ -19,7 +20,7 @@
 #' Jenness, J.S., 2004. Calculating landscape surface area from digital elevation models. Wildlife Society Bulletin 32, 829-839.
 #' @export
 
-SurfaceArea<- function(r, filename=NULL, overwrite=FALSE, wopt=list()){
+SurfaceArea<- function(r, na.rm=FALSE, filename=NULL, overwrite=FALSE, wopt=list()){
   og_class<- class(r)[1]
   if(og_class=="RasterLayer"){
     r<- terra::rast(r) #Convert to SpatRaster
@@ -37,7 +38,7 @@ SurfaceArea<- function(r, filename=NULL, overwrite=FALSE, wopt=list()){
   if(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){
     warning("Coordinate system may be Lat/Lon. Please ensure that the coordinate system is projected with elevation/depth units matching map units.")
   }
-  SA<- terra::focalCpp(r, w=c(3,3), fun = C_SurfaceArea,  x_res = terra::res(r)[1], y_res = terra::res(r)[2], fillvalue=NA, wopt=wopt)
+  SA<- terra::focalCpp(r, w=c(3,3), fun = C_SurfaceArea,  x_res = terra::res(r)[1], y_res = terra::res(r)[2], na_rm=na.rm, fillvalue=NA, wopt=wopt)
   names(SA)<- "SA"
   #Return
   if(og_class =="RasterLayer"){