Minor improvements

- use terra::mask to improve memory safety
- Putis.lonlat within isTRUE due to update in terra that returns NA instead of FALSE if no CRS
- add more detail to ReadME
- add documentation for curvature helper functions

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: MultiscaleDEM
 Title: Calculates multi-scale geomorphometric terrain attributes from regularly gridded DEM/bathymetry rasters
-Version: 0.2
+Version: 0.2.1
 Authors@R: 
      c(
     person(given = "Alexander",

R/AdjSD.R

@@ -3,7 +3,7 @@
 #' Calculates standard deviation of bathymetry (a measure of rugosity). Using a sliding rectangular window a plane is fit to the data and the standard deviation of the residuals is calculated.
 #' @param r DEM as a SpatRaster or RasterLayer in a projected coordinate system where map units match elevation/depth units
 #' @param w A vector of length 2 specifying the dimensions of the rectangular window to use where the first number is the number of rows and the second number is the number of columns. Window size must be an odd number.
-#' @param na.rm A logical vector indicating whether or not to remove NA values before calculations
+#' @param na.rm A logical indicating whether or not to remove NA values before calculations
 #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
 #' @return a SpatRaster or RasterLayer of adjusted rugosity
 #' @import terra
@@ -23,7 +23,7 @@ AdjSD<- function(r, w=c(3,3), na.rm=FALSE, include_scale=FALSE){
   if(terra::nlyr(r)!=1){
     stop("Error: Input raster must be one layer.")
   }
-  if(terra::is.lonlat(r, perhaps=FALSE)){
+  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(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){

R/RDMV.R

@@ -3,7 +3,7 @@
 #' Calculates Relative Difference from Mean Value (RDMV). RDMV = (focal_value - local_mean)/local_range
 #' @param r DEM as a SpatRaster, RasterLayer, RasterStack, or RasterBrick
 #' @param w A vector of length 2 specifying the dimensions of the rectangular window to use where the first number is the number of rows and the second number is the number of columns. Window size must be an odd number. Default is 3x3.
-#' @param na.rm A logical vector indicating whether or not to remove NA values before calculations
+#' @param na.rm A logical indicating whether or not to remove NA values before calculations
 #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
 #' @return a SpatRaster or RasterLayer
 #' @import terra

R/SAPA.R

@@ -32,7 +32,7 @@ SAPA<- function(r, w = 1, slope_correction=TRUE, include_scale=FALSE, slope_laye
   if(terra::nlyr(r)!=1){
     stop("Error: Input raster must be one layer.")
   }
-  if(terra::is.lonlat(r, perhaps=FALSE)){
+  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(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){

R/SlpAsp.R

@@ -53,7 +53,7 @@ SlpAsp <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c("slop
   if(terra::nlyr(r)!=1){
     stop("Error: Input raster must be one layer.")
   }
-  if(terra::is.lonlat(r, perhaps=FALSE)){
+  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(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){
@@ -171,43 +171,37 @@ SlpAsp <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c("slop
 
   if("slope" %in% needed_metrics){
     slope.k<- (atan(sqrt((dz.dx^2)+(dz.dy^2))))
-    if(mask_aspect){
-      slp0_idx<- slope.k==0
-      }
     if(unit=="degrees" & ("slope" %in% metrics)){
       slope.k<- slope.k * (180/pi)
       }
-
     names(slope.k)<- "slope"
     out<- c(out, slope.k, warn=FALSE)
   }
 
   if("aspect" %in% needed_metrics){
     aspect.k<- terra::app(atan2(dz.dy, -dz.dx), fun = convert_aspect) #aspect relative to North
-    # aspect.k<- (pi/2) - atan2(dz.dy, -dz.dx) #aspect relative to North
-    # aspect.k[aspect.k < 0]<- aspect.k[aspect.k < 0] + (2*pi) #Constrain in range of 0 - 2 pi
 
     if("eastness" %in% needed_metrics){
       eastness.k<- sin(aspect.k)
       if(mask_aspect){
-        eastness.k[slp0_idx]<- 0 #Set eastness to 0 where slope is zero
-      }
+        eastness.k<- terra::mask(eastness.k, mask= slope.k, maskvalues = 0, updatevalue = 0) #Set eastness to 0 where slope is zero
+        }
       names(eastness.k)<- "eastness"
       out<- c(out, eastness.k, warn=FALSE)
     }
 
     if("northness" %in% needed_metrics){
       northness.k<- cos(aspect.k)
       if(mask_aspect){
-        northness.k[slp0_idx]<- 0 #Set northenss to 0 where slope is zero
-      }
+        northness.k<- terra::mask(northness.k, mask= slope.k, maskvalues = 0, updatevalue = 0) #Set northenss to 0 where slope is zero
+        }
       names(northness.k)<- "northness"
       out<- c(out, northness.k, warn=FALSE)
     }
 
     if(mask_aspect){
-      aspect.k[slp0_idx]<- NA_real_ #Set aspect to undefined where slope is zero
-    }
+      aspect.k<- terra::mask(aspect.k, mask= slope.k, maskvalues = 0, updatevalue = NA) #Set aspect to undefined where slope is zero
+      }
     if(unit=="degrees"){
       aspect.k<- aspect.k * (180/pi)
       }

R/SurfaceArea.R

@@ -22,7 +22,7 @@ SurfaceArea<- function(r, expand=FALSE){
   if(terra::nlyr(r)!=1){
     stop("Error: Input raster must be one layer.")
   }
-  if(terra::is.lonlat(r, perhaps=FALSE)){
+  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(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){

R/TPI.R

@@ -3,7 +3,7 @@
 #' Calculates Topographic Position Index (TPI). This is the value of the focal pixel minus the mean of the surrounding pixels (i.e. local mean but excluding the value of the focal pixel).
 #' @param r DEM as a SpatRaster or RasterLayer
 #' @param w A vector of length 2 specifying the dimensions of the rectangular window to use where the first number is the number of rows and the second number is the number of columns. Window size must be an odd number. Default is 3x3.
-#' @param na.rm A logical vector indicating whether or not to remove NA values before calculations
+#' @param na.rm A logical indicating whether or not to remove NA values before calculations
 #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
 #' @return a SpatRaster or RasterLayer
 #' @import terra

R/VRM.R

@@ -3,7 +3,7 @@
 #' Implementation of the Sappington et al., (2007) vector ruggedness measure, modified from Evans (2021). 
 #' @param r DEM as a SpatRaster or RasterLayer
 #' @param w A vector of length 2 specifying the dimensions of the rectangular window to use where the first number is the number of rows and the second number is the number of columns. Window size must be an odd number.
-#' @param na.rm A logical vector indicating whether or not to remove NA values before calculations
+#' @param na.rm A logical indicating whether or not to remove NA values before calculations
 #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
 #' @return a RasterLayer
 #' @import terra

R/WoodEvans.R

@@ -52,7 +52,7 @@ convert_aspect2<- function(aspect){
 #' @param metrics Character vector specifying which terrain attributes to return. The default is to return all available metrics, c("qslope", "qaspect", "qeastness", "qnorthness", "profc", "planc", "meanc", "maxc", "minc", "longc", "crosc", "features"). Slope, aspect, eastness, and northness are preceded with a 'q' to differentiate them from the measures calculated by SlpAsp() where the 'q' indicates that a quadratic surface was used for the calculation. 'profc' is the profile curvature, 'planc' is the plan curvature, 'meanc' is the mean curvature, 'minc' is minimum curvature, 'longc' is longitudinal curvature, crosc is cross-sectional curvature, and 'features' are morphometric features. See details.
 #' @param slope_tolerance Slope tolerance that defines a 'flat' surface (degrees; default = 1.0). Relevant for the features layer.
 #' @param curvature_tolerance Curvature tolerance that defines 'planar' surface (default = 0.0001). Relevant for the features layer.
-#' @param na.rm Logical vector indicating whether or not to remove NA values before calculations.
+#' @param na.rm Logical indicating whether or not to remove NA values before calculations.
 #' @param include_scale Logical indicating whether to append window size to the layer names (default = FALSE).
 #' @param mask_aspect Logical. If TRUE (default), aspect will be set to NA and northness and eastness will be set to 0 when slope = 0. If FALSE, aspect is set to 270 degrees or 3*pi/2 radians ((-pi/2)- atan2(0,0)+2*pi) and northness and eastness will be calculated from this.
 #' @param return_params Logical indicating whether to return Wood/Evans regression parameters (default = FALSE).
@@ -86,7 +86,7 @@ WoodEvans<- function(r, w=c(3,3), unit= "degrees", metrics= c("qslope", "qaspect
   if(terra::nlyr(r)!=1){
     stop("Error: Input raster must be one layer.")
     }
-  if(terra::is.lonlat(r, perhaps=FALSE)){
+  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(terra::is.lonlat(r, perhaps=TRUE, warn=FALSE)){
@@ -117,6 +117,10 @@ WoodEvans<- function(r, w=c(3,3), unit= "degrees", metrics= c("qslope", "qaspect
     needed_metrics<- c(needed_metrics, "qaspect")
   }
 
+  if(mask_aspect & ("qaspect" %in% needed_metrics) & (!("qslope" %in% needed_metrics))){
+    needed_metrics<- c(needed_metrics, "qslope")
+  }
+
   if("features" %in% needed_metrics){
     if(!("qslope" %in% needed_metrics)){needed_metrics<- c(needed_metrics, "qslope")} 
     if(!("planc" %in% needed_metrics)){needed_metrics<- c(needed_metrics, "planc")} 
@@ -164,14 +168,13 @@ WoodEvans<- function(r, w=c(3,3), unit= "degrees", metrics= c("qslope", "qaspect
   if("qaspect" %in% needed_metrics){
     asp<- terra::app(atan2(params$e,params$d), fun = convert_aspect2) #Shift aspect so north is zero
     if (mask_aspect){
-      slp0_idx<- (params$d== 0 &  params$e== 0) #mask indicating when d ane e are 0 (slope is 0)
-      asp[slp0_idx]<- NA_real_
+      asp<- terra::mask(asp, mask= slp, maskvalues = 0, updatevalue = NA) #Set aspect to undefined where slope is zero
     }
 
     if("qeastness" %in% needed_metrics){
       eastness<- sin(asp)
       if (mask_aspect){
-        eastness[slp0_idx]<- 0
+        eastness<- terra::mask(eastness, mask= slp, maskvalues = 0, updatevalue = 0) #Set eastness to 0 where slope is zero
         }
       names(eastness)<- "qeastness"
       out<- c(out, eastness, warn=FALSE)
@@ -180,8 +183,8 @@ WoodEvans<- function(r, w=c(3,3), unit= "degrees", metrics= c("qslope", "qaspect
     if("qnorthness" %in% needed_metrics){
       northness<- cos(asp)
       if (mask_aspect){
-        northness[slp0_idx]<- 0
-      }
+        northness<- terra::mask(northness, mask= slp, maskvalues = 0, updatevalue = 0) #Set northness to 0 where slope is zero
+        }
       names(northness)<- "qnorthness"
       out<- c(out, northness, warn=FALSE)
     }
@@ -202,21 +205,21 @@ WoodEvans<- function(r, w=c(3,3), unit= "degrees", metrics= c("qslope", "qaspect
 
   if("profc" %in% needed_metrics){
     profc<- terra::lapp(params, fun = kns)
-    profc[mask_raster]<- 0
+    profc<- terra::mask(profc, mask= mask_raster, maskvalues = 1, updatevalue = 0) #Set curvature to 0 where all values are the same
     names(profc)<- "profc"
     out<- c(out, profc, warn=FALSE)
   }
 
   if("planc" %in% needed_metrics){
     planc<- terra::lapp(params, fun = knc)
-    planc[mask_raster]<- 0
+    planc<- terra::mask(planc, mask= mask_raster, maskvalues = 1, updatevalue = 0) #Set curvature to 0 where all values are the same
     names(planc)<- "planc"
     out<- c(out, planc, warn=FALSE)
   }
 
   if("twistc" %in% needed_metrics){
     twistc<- terra::lapp(params, fun = tgc)
-    twistc[mask_raster]<- 0
+    twistc<- terra::mask(twistc, mask= mask_raster, maskvalues = 1, updatevalue = 0) #Set curvature to 0 where all values are the same
     names(twistc)<- "twistc"
     out<- c(out, twistc, warn=FALSE)
     }

R/curvatures.R

@@ -1,39 +1,150 @@
+#' Calculate normal contour curvature
+#'
+#' Calculate normal contour curvature (kn)c, which is the principal representative of the plan curvature group based on regression coefficients from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
 knc<- function(a,b,c,d,e,f){
   #Planc: Normal Contour Curvature
   out<- -(2*a*(e^2) - 2*c*d*e + 2*b*d^2)/((d^2+e^2) * sqrt(1+d^2+e^2))
   return(out)
 }
 
+#' Calculate normal slope line curvature
+#'
+#' Calculate normal slope line curvature (kn)s, which is the principal representative of the profile curvature group based on regression coefficients from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
+
 kns<- function(a,b,c,d,e,f){
   #Profc: Normal Slope Line Curvature
   out<- (-2 * (a*d^2 + b*e^2 + c*d*e)) / ((e^2 + d^2)*(1 + e^2 + d^2)^1.5)
   return(out)
 }
 
+#' Calculate contour geodesic torsion
+#'
+#' Calculate contour geodesic torsion (tg)c, which is the principal representative of the twisting curvature group based on regression coefficients from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
+
 tgc<- function(a,b,c,d,e,f){
   #TwistC: Contour geodesic torsion
   out<- (d*e*(2*a-2*b) - c*(d^2-e^2))/((d^2+e^2)*(1+d^2+e^2))
   return(out)
 }
 
+#' Calculate mean curvature
+#'
+#' Calculate mean curvature, kmean, from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
+
 kmean<- function(a,b,c,d,e,f){
   #Mean Curvature
   out<- -((1+e^2)*2*a - 2*c*d*e + (1+d^2)*2*b) / (2*sqrt((1+d^2+e^2)^3))
   return(out)
 }
 
+#' Calculate unsphericity curvature
+#'
+#' Calculate  unsphericity curvature, ku, from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
+
 ku<- function(a,b,c,d,e,f){
   #unsphericity curvature
   out<- sqrt((((1+e^2)*2*a - 2*c*d*e +(1+d^2)*2*b) / (2*sqrt((1+d^2+e^2)^3)))^2 - ((2*a*2*b-c^2)/(1+d^2+e^2)^2))
   return(out)
 }
 
+#' Calculate min curvature
+#'
+#' Calculate min curvature, kmin, from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
+
 kmin<- function(a,b,c,d,e,f){
   #Min Curvature
   out<- kmean(a,b,c,d,e)-ku(a,b,c,d,e)
   return(out)
 }
 
+#' Calculate max curvature
+#'
+#' Calculate max curvature, kmax, from the equation Z =ax^2+by^2+cxy+dx+e+f.
+#' @param a regression coefficient
+#' @param b regression coefficient
+#' @param c regression coefficient
+#' @param d regression coefficient
+#' @param e regression coefficient
+#' @param f regression intercept
+#' @references
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' 
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' 
+#' Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414. https://doi.org/10.1016/j.earscirev.2020.103414
+
 kmax<- function(a,b,c,d,e,f){
   #Max Curvature
   out<- kmean(a,b,c,d,e)+ku(a,b,c,d,e)