Optimize C++ code!

DESCRIPTION

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

NAMESPACE

@@ -4,6 +4,7 @@ export(AdjSD)
 export(BPI)
 export(DMV)
 export(DirSlp)
+export(Pfit)
 export(Qfit)
 export(RIE)
 export(RelPos)

R/AdjSD.R

@@ -10,9 +10,7 @@
 #' @param wopt list with named options for writing files as in writeRaster
 #' @return a SpatRaster or RasterLayer of adjusted rugosity
 #' @examples
-#' r<- rast(volcano, extent= ext(2667400, 2667400 + 
-#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), 
-#' crs = "EPSG:27200")
+#' r<- erupt()
 #' adjsd<- AdjSD(r, w=c(5,5), na.rm = TRUE)
 #' plot(adjsd)
 #' @import terra

R/BPI.R

@@ -12,9 +12,7 @@
 #' @param wopt List with named options for writing files as in writeRaster.
 #' @return A SpatRaster or RasterLayer.
 #' @examples
-#' r<- rast(volcano, extent= ext(2667400, 2667400 + 
-#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), 
-#' crs = "EPSG:27200")
+#' r<- erupt()
 #' bpi<- BPI(r, w = c(2,4), stand= "none", unit = "cell", na.rm = TRUE)
 #' plot(bpi)
 #' @import terra

R/DMV.R

@@ -13,9 +13,7 @@
 #' @param wopt List with named options for writing files as in writeRaster.
 #' @return a SpatRaster or RasterLayer.
 #' @examples
-#' r<- rast(volcano, extent= ext(2667400, 2667400 + 
-#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), 
-#' crs = "EPSG:27200")
+#' r<- erupt()
 #' dmv<- DMV(r, w=c(5,5), shape= "rectangle", stand="range", na.rm = TRUE)
 #' plot(dmv)
 #' @import terra

R/DirSlope.R

@@ -19,9 +19,11 @@
 #' r<- erupt()
 #' dz1<- SlpAsp(r, metrics = c("dz.dx", "dz.dy"))
 #' dz2<- Qfit(r, metrics = c(), return_params = TRUE, as_derivs=TRUE)
+#' dz3<- Pfit(r, metrics = c("dz.dx", "dz.dy"))
 #' dirslp1<- DirSlp(alpha = 45, dz.dx= dz1$dz.dx, dz.dy= dz1$dz.dy)
 #' dirslp2<- DirSlp(alpha = 45, dz.dx= dz2$zx, dz.dy= dz2$zy)
-#' @details dz.dx can be calculated at a specified scale via `SlpAsp`, `Qfit` (zx and zy), or from an existing layer calculated by another program.
+#' dirslp3<- DirSlp(alpha = 45, dz.dx= dz3$dz.dx, dz.dy= dz3$dz.dy)
+#' @details dz.dx and dz.dy can be calculated at a specified scale via `SlpAsp`, `Pfit`, `Qfit` (zx and zy), or from an existing layer calculated by another program.
 #' @references
 #' Neteler, M., & Mitasova, H. (2008). Open source GIS: A GRASS GIS approach (3rd ed.). Springer.
 #' @import terra

R/Pfit.R

@@ -1,24 +1,36 @@
-#' Calculates multiscale slope and aspect using a local planar fit (IN DEVELOPMENT: DO NOT USE)
+#' Calculates multiscale slope and aspect using a local planar fit.
 #'
-#' Calculates multiscale slope and aspect using a local planar fit (IN DEVELOPMENT: DO NOT USE)
+#' Calculates multiscale slope and aspect of a DTM over a sliding rectangular window using a using a local planar fit to the surface (Sharpnack and Akin 1969).
 #' @param r DTM as a SpatRaster (terra) or RasterLayer (raster) in a projected coordinate system where map units match elevation/depth units (up is assumed to be north for calculations of aspect, northness, and eastness).
 #' @param w 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 unit "degrees" or "radians".
-#' @param na.rm Logical indicating whether or not to remove NA values before calculations.#' 
+#' @param metrics Character vector specifying which terrain attributes to return. The default is to return c("pslope", "paspect", "peastness", and "pnorthness"). These are preceded with a 'p' to differentiate them from the measures calculated by SlpAsp() and 'Qfit' where the 'p' indicates that a planar surface was used for the calculation. Additional measures available include "dz.dx" and "dz.dy" which are the x and y components of slope respectively.
+#' @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 filename character Output filename. Can be a single filename, or as many filenames as there are layers to write a file for each layer
 #' @param overwrite logical. If TRUE, filename is overwritten (default is FALSE).
 #' @param wopt list with named options for writing files as in writeRaster
 #' @return a SpatRaster (terra) or RasterStack/RasterLayer (raster)
+#' @examples
+#' r<- erupt()
+#' pmetrics<- Pfit(r, w = c(5,5), unit = "degrees", na.rm = TRUE)
+#' plot(pmetrics)
 #' @import terra
 #' @importFrom raster raster
 #' @importFrom raster stack
 #' @importFrom raster writeRaster
+#' @details
+#' If only first order derivatives are needed, Pfit is faster than Qfit and should provide equivalent results to Qfit for first order derivatives (Jones, 1998) when na.rm=FALSE and approximately the same results otherwise. 
+#' 
 #' @references
-#' Sharpnack, D. A. and Akin, G., 1969. An algorithm for computing slope and aspect from elevations. Photogrammetric Engineering 35(3), 247-248.
+#' Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.
+#' Jones, K. H. (1998). A comparison of algorithms used to compute hill slope as a property of the DEM. Computers & Geosciences, 24(4), 315–323. https://doi.org/10.1016/S0098-3004(98)00032-6
+#' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+#' @export
 
-Pfit<- function(r, w=c(3,3), unit= "degrees",na.rm=FALSE, include_scale=FALSE, filename=NULL, overwrite=FALSE, wopt=list()){
-
+Pfit<- function(r, w=c(3,3), unit= "degrees", metrics = c("pslope", "paspect", "peastness", "pnorthness"), na.rm=FALSE, include_scale=FALSE, mask_aspect=TRUE, filename=NULL, overwrite=FALSE, wopt=list()){
+  all_metrics<- c("pslope", "paspect", "peastness", "pnorthness", "dz.dx", "dz.dy")
   og_class<- class(r)[1]
   if(og_class=="RasterLayer"){
     r<- terra::rast(r) #Convert to SpatRaster
@@ -51,34 +63,141 @@ Pfit<- function(r, w=c(3,3), unit= "degrees",na.rm=FALSE, include_scale=FALSE, f
   if (!any(unit==c("degrees", "radians"))){
     stop("unit must be 'degrees' or 'radians'")
   }
+
+  # Increase tolerance of inputs to metrics
+  metrics<- tolower(metrics) #Make all lowercase
+  metrics[metrics=="slope"]<- "pslope" #replace slope with pslope
+  metrics[metrics=="aspect"]<- "paspect" #replace aspect with paspect
+  metrics[metrics=="eastness"]<- "peastness" #replace eastness with peastness
+  metrics[metrics=="northness"]<- "pnorthness" #replace northness with pnorthness
 
+  if (any(!(metrics %in% all_metrics))){
+    stop("Error: Invlaid metric. Valid metrics include 'pslope', 'paspect', 'peastness', 'pnorthness', 'dz.dx', and 'dz.dy'.")
+  }
+
+  needed_metrics<- metrics
+
+  if(any(c("peastness", "pnorthness") %in% needed_metrics) & (!("paspect" %in% needed_metrics))){
+    needed_metrics<- c(needed_metrics, "paspect")
+  }
+
+  if(mask_aspect & ("paspect" %in% needed_metrics) & (!("pslope" %in% needed_metrics))){
+    needed_metrics<- c(needed_metrics, "pslope")
+  }
+
+  if(any(c("pslope", "paspect", "peastness", "pnorthness") %in% needed_metrics) & !("dz.dx" %in% needed_metrics)){
+    needed_metrics<- c(needed_metrics, "dz.dx")
+  }
+
+  if(any(c("pslope", "paspect", "peastness", "pnorthness") %in% needed_metrics) & !("dz.dy" %in% needed_metrics)){
+    needed_metrics<- c(needed_metrics, "dz.dy")
+  }
+
   #Define local coordinate system of window
 
   x_mat<- matrix(data = seq(from = (-xres(r) * floor(w[2]/2)), to = (xres(r) * floor(w[2]/2)), length.out = w[2]), nrow = w[1], ncol=w[2], byrow=TRUE)
-  x_mat[-c(1, nrow(x_mat)), -c(1, ncol(x_mat))] <- NA
   x<- as.vector(t(x_mat)) #Transpose for focal
 
   y_mat<- matrix(data = seq(from = (yres(r) * floor(w[1]/2)), to = (-yres(r) * floor(w[1]/2)), length.out = w[1]), nrow = w[1], ncol=w[2])
-  y_mat[-c(1, nrow(y_mat)), -c(1, ncol(y_mat))] <- NA
   y<- as.vector(t(y_mat)) #Transpose for focal
 
-  idx<- !is.na(x)
+  # forcing through center doesn't affect slope (Jones, 1998) so don't bother estimating the extra coefficient unless na.rm=TRUE
+  # If na.rm = TRUE, you want to be able to get the slope even if the central cell is NA
+  X<- cbind(x, y) # Z = dx+ey (no intercept)
 
-  X<- cbind(x, y, 1) #Z = dx+ey+f
-  X<- X[idx,]
+  if(na.rm){
+    X<- cbind(X,1)
+  } # Z = dx+ey+f (with intercept)
 
   if(!na.rm){
     Xt<- t(X)
     XtX_inv<- solve(Xt %*% X)
   }
 
   # Calculate Regression Parameters
-
-  params<- terra::focalCpp(r, w=w, fun = C_Pfit1_narmF, X= X, Xt= Xt, XtX_inv= XtX_inv, idx=idx, fillvalue=NA, wopt=wopt)
-  slp<- terra::math(terra::math(params$d^2 + params$e^2, fun="sqrt", wopt=wopt), fun="atan", wopt=wopt)
-  if(unit=="degrees"){
-    slp<- slp*(180/pi)
+  if(na.rm){
+    params<- terra::focalCpp(r, w=w, fun = C_Qfit1_narmT, X_full= X, return_intercept = FALSE, fillvalue=NA, wopt=wopt)
+  } else{
+    params<- terra::focalCpp(r, w=w, fun = C_Qfit2_narmF, X= X, Xt= Xt, XtX_inv= XtX_inv, fillvalue=NA, wopt=wopt)
+  }
+  names(params)<- c("dz.dx", "dz.dy")
+
+  out<- terra::rast() #Initialize output
+
+  if("dz.dx" %in% needed_metrics){
+    out<- c(out, params$dz.dx, warn=FALSE)
+  }
+
+  if("dz.dy" %in% needed_metrics){
+    out<- c(out, params$dz.dy, warn=FALSE)
+  }
+
+  #Use regression parameters to calculate slope and aspect
+  if("pslope" %in% needed_metrics){
+    slp<- terra::math(terra::math(params$dz.dx^2 + params$dz.dy^2, fun="sqrt", wopt=wopt), fun="atan", wopt=wopt)
+    if(unit=="degrees"){
+      slp<- slp*(180/pi)
+    }
+    names(slp)<- "pslope"
+    out<- c(out, slp, warn=FALSE)
+  }
+
+  if("paspect" %in% needed_metrics){
+    asp<- (-pi/2) - terra::atan_2(params$dz.dy,params$dz.dx, wopt=wopt) # aspect relative to North
+    asp<- ifel(asp < 0, yes = asp+(2*pi), no= asp, wopt=wopt) # Constrain range so between 0 and 2pi
+
+    if (mask_aspect & ("paspect" %in% metrics)){
+      asp<- terra::mask(asp, mask= slp, maskvalues = 0, updatevalue = NA, wopt=wopt) #Set aspect to undefined where slope is zero (doesn't need to be needed_metrics since also mask northness and eastness as 0).
+    }
+
+    if("peastness" %in% needed_metrics){
+      eastness<- sin(asp)
+      if (mask_aspect){
+        eastness<- terra::mask(eastness, mask= slp, maskvalues = 0, updatevalue = 0, wopt=wopt) #Set eastness to 0 where slope is zero
+      }
+      names(eastness)<- "peastness"
+      out<- c(out, eastness, warn=FALSE)
+    }
+
+    if("pnorthness" %in% needed_metrics){
+      northness<- cos(asp)
+      if (mask_aspect){
+        northness<- terra::mask(northness, mask= slp, maskvalues = 0, updatevalue = 0, wopt=wopt) #Set northness to 0 where slope is zero
+      }
+      names(northness)<- "pnorthness"
+      out<- c(out, northness, warn=FALSE)
+    }
+    if(unit=="degrees"){
+      asp<- asp*180/pi
+    }
+    names(asp)<- "paspect"
+    out<- c(out, asp, warn=FALSE)
+  }
+
+  out<- terra::subset(out, metrics, wopt=wopt) #Subset needed metrics to requested metrics in proper order
+
+  if(include_scale){names(out)<- paste0(names(out), "_", w[1],"x", w[2])} #Add scale to layer names
+
+  #Return
+  if(og_class=="RasterLayer"){
+    if(terra::nlyr(out) > 1){
+      out<- raster::stack(out) #Convert to RasterStack
+      if(!is.null(filename)){
+        if(length(filename)==1){
+          return(raster::writeRaster(out, filename=filename, overwrite=overwrite, bylayer=FALSE))
+        } else{
+          return(raster::writeRaster(out, filename=filename, overwrite=overwrite, bylayer=TRUE))
+        }
+      }
+    } else{
+      out<- raster::raster(out)
+      if(!is.null(filename)){
+        return(raster::writeRaster(out, filename=filename, overwrite=overwrite))
+      }
+    }
+  }
+  if(!is.null(filename)){
+    return(terra::writeRaster(out, filename=filename, overwrite=overwrite, wopt=wopt))
   }
-  # asp<- #code goes here
-  return(slp)
+  return(out)
 }

R/Qfit.R

@@ -67,8 +67,8 @@ outlier_filter<- function(params, outlier_quantile, wopt=list()){
 #' @param curvature_tolerance Curvature tolerance that defines 'planar' surface (default = 0.0001). Relevant for the features layer.
 #' @param outlier_quantile A numeric vector of length two or three. If two numbers are used it specifies the lower (Q1) and upper (Q2) quantiles used for determining the interquantile range (IQR). These values should be between 0 and 1 with Q2 > Q1. An optional third number can be used to specify a the size of a regular sample to be taken which can be useful if the full dataset is too large to fit in memory. Values are considered outliers and replaced with NA if they are less than Q1-(100*IQR) or greater than Q2+(100*IQR), where IQR=Q2-Q1. The outlier filter is performed on the results of the regression parameters ('a'-'e' and 'elev') prior to calculation of subsequent terrain attributes. Note that c(0,1) will skip the outlier filtering step and can speed up computations. The default is c(0.01,0.99).
 #' @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 force_center Logical specifying whether the constrain the model through the central cell of the focal window
+#' @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).
 #' @param as_derivs Logical indicating whether parameters should be formatted as partial derivatives instead of regression coefficients (default = FALSE) (Minár et al., 2020).
@@ -77,9 +77,7 @@ outlier_filter<- function(params, outlier_quantile, wopt=list()){
 #' @param wopt list with named options for writing files as in writeRaster
 #' @return a SpatRaster (terra) or RasterStack/RasterLayer (raster)
 #' @examples
-#' r<- rast(volcano, extent= ext(2667400, 2667400 + 
-#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), 
-#' crs = "EPSG:27200")
+#' r<- erupt()
 #' qmetrics<- Qfit(r, w = c(5,5), unit = "degrees", na.rm = TRUE)
 #' plot(qmetrics)
 #' 
@@ -155,7 +153,7 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
 
   if(force_center & ("elev" %in% metrics)){
     metrics<- metrics[metrics!="elev"]
-    warning("Warning: dropping 'elev' from metrics since force_central is TRUE")
+    warning("Warning: dropping 'elev' from metrics since force_center is TRUE")
     }
   needed_metrics<- metrics
 
@@ -192,26 +190,30 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
     }
 
   # Calculate Regression Parameters
+  return_intercept<- "elev" %in% needed_metrics
+
   if((!na.rm) & (!force_center)){
-    params<- terra::focalCpp(r, w=w, fun = C_Qfit1_narmF, X= X, Xt= Xt, XtX_inv= XtX_inv, fillvalue=NA, wopt=wopt)
+    params<- terra::focalCpp(r, w=w, fun = C_Qfit1_narmF, X= X, Xt= Xt, XtX_inv= XtX_inv, return_intercept = return_intercept, fillvalue=NA, wopt=wopt)
     } else if(na.rm & (!force_center)){
-      params<- terra::focalCpp(r, w=w, fun = C_Qfit1_narmT, X_full= X, fillvalue=NA, wopt=wopt)
+      params<- terra::focalCpp(r, w=w, fun = C_Qfit1_narmT, X_full= X, return_intercept = return_intercept, fillvalue=NA, wopt=wopt)
     } else if((!na.rm) & force_center){
       params<- terra::focalCpp(r, w=w, fun = C_Qfit2_narmF, X= X, Xt= Xt, XtX_inv= XtX_inv, fillvalue=NA, wopt=wopt)} else{
         params<- terra::focalCpp(r, w=w, fun = C_Qfit2_narmT, X_full= X, fillvalue=NA, wopt=wopt)
       }
 
+  if(return_intercept){
+    names(params)<- c("a","b","c","d","e","elev")
+    elev<- params$elev
+    params<- params[[-6]] #drop intercept
+  } else{
+    names(params)<- c("a","b","c","d","e")
+  }
+
   # Filter outliers
   if(!all(outlier_quantile[c(1,2)]==c(0,1))){
     params <- outlier_filter(params, outlier_quantile, wopt=wopt)
     }
 
-  if(!force_center){
-    elev<- params$f
-    names(elev)<- "elev"
-    params<- params[[-6]] #drop intercept
-  }
-
   mask_raster<- terra::app(terra::math(params, fun= "abs", wopt=wopt), fun= "sum", na.rm=FALSE, wopt=wopt) == 0 # Mask of when all params are 0
 
 
@@ -236,8 +238,8 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
     asp<- (-pi/2) - terra::atan_2(params$e,params$d, wopt=wopt) # aspect relative to North
     asp<- ifel(asp < 0, yes = asp+(2*pi), no= asp, wopt=wopt) # Constrain range so between 0 and 2pi
 
-    if (mask_aspect){
-      asp<- terra::mask(asp, mask= slp, maskvalues = 0, updatevalue = NA, wopt=wopt) #Set aspect to undefined where slope is zero
+    if (mask_aspect & ("qaspect" %in% metrics)){
+      asp<- terra::mask(asp, mask= slp, maskvalues = 0, updatevalue = NA, wopt=wopt) #Set aspect to undefined where slope is zero (doesn't need to be needed_metrics since also mask northness and eastness as 0).
     }
 
     if("qeastness" %in% needed_metrics){

R/RIE.R

@@ -10,9 +10,7 @@
 #' @param wopt list with named options for writing files as in writeRaster
 #' @return a SpatRaster or RasterLayer
 #' @examples 
-#' r<- rast(volcano, extent= ext(2667400, 2667400 + 
-#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), 
-#' crs = "EPSG:27200")
+#' r<- erupt()
 #' rie<- RIE(r, w=c(5,5), na.rm = TRUE)
 #' plot(rie)
 #' @import terra