First Full Release

- Remove pre release warning
- Add sloping category to features
- wopt changes

DESCRIPTION

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

R/Qfit.R

@@ -9,23 +9,23 @@
 #' Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
 
 classify_features_ff<- function(slope_tolerance=1, curvature_tolerance=0.0001){
-  #1=PLANAR (If add slope redefine to FLAT?)
-  #2=PIT
-  #3=CHANNEL
-  #4=PASS
-  #5=RIDGE
-  #6=PEAK
-  #7=SLOPE??? (add a new category?)
+  #1=PLANAR_FLAT
+  #2=PLANAR_SLOPE
+  #3=PIT
+  #4=CHANNEL
+  #5=PASS
+  #6=RIDGE
+  #7=PEAK
   out_fun<- function(slope, planc, maxc, minc){
     dplyr::case_when(is.na(slope) ~ NA_real_,
-                     (slope > slope_tolerance) & (planc > curvature_tolerance) ~ 5,
-                     (slope > slope_tolerance) & (planc < -curvature_tolerance) ~ 3,
-                     slope > slope_tolerance ~ 1, #Maybe make this 7 (slope, new category)
-                     (maxc > curvature_tolerance) & (minc > curvature_tolerance) ~ 6,
-                     (maxc > curvature_tolerance) & (minc < -curvature_tolerance) ~ 4,
-                     maxc > curvature_tolerance ~ 5,
-                     (minc < -curvature_tolerance) & (maxc < -curvature_tolerance) ~ 2,
-                     minc < -curvature_tolerance ~ 3,
+                     (slope > slope_tolerance) & (planc > curvature_tolerance) ~ 6,
+                     (slope > slope_tolerance) & (planc < -curvature_tolerance) ~ 4,
+                     slope > slope_tolerance ~ 2, 
+                     (maxc > curvature_tolerance) & (minc > curvature_tolerance) ~ 7,
+                     (maxc > curvature_tolerance) & (minc < -curvature_tolerance) ~ 5,
+                     maxc > curvature_tolerance ~ 6,
+                     (minc < -curvature_tolerance) & (maxc < -curvature_tolerance) ~ 3,
+                     minc < -curvature_tolerance ~ 4,
                      TRUE ~ 1)}
   return(out_fun)
 }
@@ -47,8 +47,9 @@ convert_aspect2<- function(aspect){
 #'
 #' @param params regression parameters for fitted surface
 #' @param outlier_quantile vector of length 2 specifying the quantiles used for filtering outliers
+#' @param wopt list with named options for writing files as in writeRaster
 
-outlier_filter<- function(params, outlier_quantile){ 
+outlier_filter<- function(params, outlier_quantile, wopt=list()){ 
   quant <- terra::global(params, fun = quantile, probs = c(0, outlier_quantile[1], outlier_quantile[2], 1), na.rm = TRUE)
   iqr <- quant[, 3] - quant[, 2]
   outliers <- row.names(quant)[which(quant[, 1] < (quant[, 2] - 100 * iqr) | quant[, 4] > (quant[, 3] + 100 * iqr))]
@@ -60,8 +61,8 @@ outlier_filter<- function(params, outlier_quantile){
     for (i in 1:nlyr(params)) {
       outlier_mask[[i]]<- ((params[[i]] >= iq_lims[i,1]) & (params[[i]] <= iq_lims[i,2])) #0 indicates an outlier
     }
-    outlier_mask<- prod(outlier_mask) #product of zero indicates an outlier location
-    params<- terra::mask(params, mask = outlier_mask, maskvalues=0, updatevalue=NA)
+    outlier_mask<- terra::prod(outlier_mask, wopt=wopt) #product of zero indicates an outlier location
+    params<- terra::mask(params, mask = outlier_mask, maskvalues=0, updatevalue=NA, wopt=wopt)
     warning("Outliers filtered")
   }
   return(params)
@@ -96,7 +97,7 @@ outlier_filter<- function(params, outlier_quantile){
 #' To get only the regression coefficients, set "metrics=c()" and "return_params=TRUE"
 #' reg_coefs<- Qfit(r, w = c(5,5), metrics=c(), unit = "degrees", na.rm = TRUE, return_params=TRUE)
 #' plot(reg_coefs)
-#' @details This function calculates slope, aspect, eastness, northness, profile curvature, plan curvature, mean curvature, twisting curvature, maximum curvature, minimum curvature, morphometric features, and a smoothed version of the elevation surface using a quadratic surface fit from Z = aX^2+bY^2+cXY+dX+eY+f, where Z is the elevation or depth values, X and Y are the xy coordinates relative to the central cell in the focal window, and a-f are parameters to be estimated (Evans, 1980; Minár et al. 2020; Wood, 1996). For aspect, 0 degrees represents north (or if rotated, the direction that increases as you go up rows in your data) and increases clockwise. For calculations of northness (cos(asp)) and eastness (sin(asp)), up in the y direction is assumed to be north, and if this is not true for your data (e.g. you are using a rotated coordinate system), you must adjust accordingly. All curvature formulas are adapted from Minár et al 2020. Therefore all curvatures are reported in units of 1/length and have are reported according to a geographic sign convention where convex is positive and concave is negative (i.e., hills are considered convex with positive curvature values). Naming convention for curvatures is not consistent across the literature, however Minár et al (2020) has suggested a framework in which the reported measures of curvature translate to profile curvature = (kn)s, plan curvature = (kn)c, twisting curvature (τg)c, mean curvature = kmean, maximum curvature = kmax, minimum curvature = kmin. For morphometric features cross-sectional curvature (zcc) was replaced by planc (kn)c, z''min was replaced by kmax, and z''max was replaced by kmin as these are more robust ways to measures the same types of curvature (Minár et al., 2020).
+#' @details This function calculates slope, aspect, eastness, northness, profile curvature, plan curvature, mean curvature, twisting curvature, maximum curvature, minimum curvature, morphometric features, and a smoothed version of the elevation surface using a quadratic surface fit from Z = aX^2+bY^2+cXY+dX+eY+f, where Z is the elevation or depth values, X and Y are the xy coordinates relative to the central cell in the focal window, and a-f are parameters to be estimated (Evans, 1980; Minár et al. 2020; Wood, 1996). For aspect, 0 degrees represents north (or if rotated, the direction that increases as you go up rows in your data) and increases clockwise. For calculations of northness (cos(asp)) and eastness (sin(asp)), up in the y direction is assumed to be north, and if this is not true for your data (e.g. you are using a rotated coordinate system), you must adjust accordingly. All curvature formulas are adapted from Minár et al 2020. Therefore all curvatures are reported in units of 1/length and have are reported according to a geographic sign convention where convex is positive and concave is negative (i.e., hills are considered convex with positive curvature values). Naming convention for curvatures is not consistent across the literature, however Minár et al (2020) has suggested a framework in which the reported measures of curvature translate to profile curvature = (kn)s, plan curvature = (kn)c, twisting curvature (τg)c, mean curvature = kmean, maximum curvature = kmax, minimum curvature = kmin. For morphometric features cross-sectional curvature (zcc) was replaced by planc (kn)c, z''min was replaced by kmax, and z''max was replaced by kmin as these are more robust ways to measures the same types of curvature (Minár et al., 2020). Additionally, the planar feature from Wood (1996) was split into planar flat and slope depending on whether the slope threshold is exceeded or not.
 #' @import terra
 #' @importFrom raster raster
 #' @importFrom raster stack
@@ -191,7 +192,7 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
   if(force_center){
     params<- terra::focalCpp(r, w=w, fun = C_Qfit2, X_full= X, na_rm=na.rm, fillvalue=NA, wopt=wopt)
     if(!all(outlier_quantile==c(0,1))){
-      params <- outlier_filter(params, outlier_quantile)
+      params <- outlier_filter(params, outlier_quantile, wopt=wopt)
       }
     } else{
     params<- terra::focalCpp(r, w=w, fun = C_Qfit1, X_full= X, na_rm=na.rm, fillvalue=NA, wopt=wopt)
@@ -211,7 +212,7 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
 
   #Use regression parameters to calculate slope and aspect
   if("qslope" %in% needed_metrics){
-    slp<- atan(sqrt(params$d^2 + params$e^2))
+    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)
     } else{
@@ -222,7 +223,7 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
   }
 
   if("qaspect" %in% needed_metrics){
-    asp<- terra::app(atan2(params$e,params$d), fun = convert_aspect2, wopt=wopt) #Shift aspect so north is zero
+    asp<- terra::app(terra::atan_2(params$e,params$d, wopt=wopt), fun = convert_aspect2, wopt=wopt) #Shift aspect so north is zero
     if (mask_aspect){
       asp<- terra::mask(asp, mask= slp, maskvalues = 0, updatevalue = NA, wopt=wopt) #Set aspect to undefined where slope is zero
     }
@@ -302,7 +303,7 @@ Qfit<- function(r, w=c(3,3), unit= "degrees", metrics= c("elev", "qslope", "qasp
   if("features" %in% needed_metrics){
     classify_features<- classify_features_ff(slope_tolerance, curvature_tolerance) #Define classification function based on slope and curvature tolerance
     features<- terra::lapp(c(slp, planc, maxc, minc), fun = classify_features, wopt=wopt)
-    levels(features)<- data.frame(ID=1:6, features = c("Planar", "Pit", "Channel", "Pass", "Ridge", "Peak"))
+    levels(features)<- data.frame(ID=1:7, features = c("Planar_Flat","Planar_Slope", "Pit", "Channel", "Pass", "Ridge", "Peak"))
     names(features)<- "features"
     out<- c(out, features, warn=FALSE)
   }

R/SAPA.R

@@ -57,8 +57,8 @@ SAPA<- function(r, w = 1, slope_correction=TRUE, include_scale=FALSE, slope_laye
     is_native<- TRUE} else{
       is_native<- FALSE} #Indicate whether SAPA is calculated at native scale
 
-  x_res<- terra::res(r)[1]
-  y_res<- terra::res(r)[2]
+  x_res<- terra::xres(r)
+  y_res<- terra::yres(r)
 
   SA<- SurfaceArea(r, wopt=wopt)
 
@@ -74,7 +74,7 @@ SAPA<- function(r, w = 1, slope_correction=TRUE, include_scale=FALSE, slope_laye
       }
     }
   PA<- (x_res*w[2]) * (y_res*w[1])
-  if(slope_correction){PA<- PA/cos(slope_layer)}#Planar area corrected for slope
+  if(slope_correction){PA<- PA/terra::math(slope_layer, fun="cos", wopt=wopt)}#Planar area corrected for slope
 
   sapa<- SA/PA
   names(sapa)<- "sapa"

R/SlpAsp.R

@@ -175,14 +175,14 @@ SlpAsp <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c("slop
     dz.dy.b <- terra::focal(r, yb.mat, fun=sum, na.rm=FALSE, wopt=wopt)
 
     #calculate dz/dx and dz/dy using the components. 2*j is the run: (2 sides)*(length of each side)
-    dz.dx <- (dz.dx.r-dz.dx.l)/(2*jx*terra::res(r)[1])
-    dz.dy <- (dz.dy.b-dz.dy.t)/(2*jy*terra::res(r)[2])
+    dz.dx <- (dz.dx.r-dz.dx.l)/(2*jx*terra::xres(r))
+    dz.dy <- (dz.dy.b-dz.dy.t)/(2*jy*terra::yres(r))
   }
 
   out<- terra::rast() #initialize output
 
   if("slope" %in% needed_metrics){
-    slope.k<- (atan(sqrt((dz.dx^2)+(dz.dy^2))))
+    slope.k<- terra::math(terra::math(dz.dx^2 + dz.dy^2, fun= "sqrt", wopt=wopt), fun= "atan", wopt=wopt)
     if(unit=="degrees" & ("slope" %in% metrics)){
       slope.k<- slope.k * (180/pi)
       }
@@ -191,10 +191,10 @@ SlpAsp <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c("slop
   }
 
   if("aspect" %in% needed_metrics){
-    aspect.k<- terra::app(atan2(dz.dy, -dz.dx), fun = convert_aspect, wopt=wopt) #aspect relative to North
+    aspect.k<- terra::app(terra::atan_2(dz.dy, -dz.dx, wopt=wopt), fun = convert_aspect, wopt=wopt) #aspect relative to North
 
     if("eastness" %in% needed_metrics){
-      eastness.k<- sin(aspect.k)
+      eastness.k<- terra::math(aspect.k, fun="sin", wopt=wopt)
       if(mask_aspect){
         eastness.k<- terra::mask(eastness.k, mask= slope.k, maskvalues = 0, updatevalue = 0, wopt=wopt) #Set eastness to 0 where slope is zero
         }
@@ -203,7 +203,7 @@ SlpAsp <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c("slop
     }
 
     if("northness" %in% needed_metrics){
-      northness.k<- cos(aspect.k)
+      northness.k<- terra::math(aspect.k, fun="cos", wopt=wopt)
       if(mask_aspect){
         northness.k<- terra::mask(northness.k, mask= slope.k, maskvalues = 0, updatevalue = 0, wopt=wopt) #Set northenss to 0 where slope is zero
         }

R/VRM.R

@@ -45,10 +45,10 @@ VRM<- function(r, w, na.rm = FALSE, include_scale=FALSE, filename=NULL, overwrit
     stop("Error: w must be greater or equal to 3 in at least one dimension")
   }
   sa <- terra::terrain(r, v=c("slope", "aspect"), unit="radians", neighbors=8, wopt=wopt) 					
-  sin.slp <- sin(sa$slope) # xyRaster 
-  cos.slp <- cos(sa$slope) # zRaster 
-  sin.asp <- sin(sa$aspect) * sin.slp # yRaster
-  cos.asp <- cos(sa$aspect) * sin.slp # xRaster  
+  sin.slp <- terra::math(sa$slope, fun="sin", wopt=wopt) # xyRaster 
+  cos.slp <- terra::math(sa$slope, fun="cos", wopt=wopt) # zRaster 
+  sin.asp <- terra::math(sa$aspect, fun="sin", wopt=wopt) * sin.slp # yRaster
+  cos.asp <- terra::math(sa$aspect, fun="cos", wopt=wopt) * sin.slp # xRaster  
   x.sum <- terra::focal(sin.asp, w = w, fun=sum, na.rm=na.rm, wopt=wopt)
   y.sum <- terra::focal(cos.asp, w = w, fun=sum, na.rm=na.rm, wopt=wopt) 
   z.sum <- terra::focal(cos.slp, w = w, fun=sum, na.rm=na.rm, wopt=wopt)

README.Rmd

@@ -8,8 +8,6 @@ output:
 ---
 # MultiscaleDEM
 
-## THIS IS A PRE-RELEASE. IT IS STILL BEING TESTED AND FUNCTIONS MAY RECEIVE MAJOR CHANGES.
-
 [![DOI](https://zenodo.org/badge/353158828.svg)](https://zenodo.org/badge/latestdoi/353158828)
 
 
@@ -29,9 +27,9 @@ This package calculates multi-scale geomorphometric terrain attributes from regu
 ## Install and Load Package
 If you don't already have remotes installed, use the code `install.packages("remotes")`
 
-Then to install this package use the code `remotes::install_github("ailich/MultiscaleDEM")` (you may need to install Rtools using the instructions found [here](https://cran.r-project.org/bin/windows/Rtools/)).
+Then to install this package use the code `remotes::install_github("ailich/MultiscaleDEM")` (If you are using Windows, you may need to install Rtools using the instructions found [here](https://cran.r-project.org/bin/windows/Rtools/)).
 
-This package relies on the `terra` package for handling of spatial raster data. To install the development version of `terra`, use `install.packages('terra', repos='https://rspatial.r-universe.dev')`. This package is also backwards compatible with the `raster` package. To install the development version of `raster` use `install.packages('raster', repos='https://rspatial.r-universe.dev')`
+This package relies on the `terra` package for handling of spatial raster data. The CRAN version of `terra` may not have all necessary functions, so it is recommended to install the development version of `terra`. The easiest way to install the development version of `terra` on Windows or Mac is to use `install.packages('terra', repos='https://rspatial.r-universe.dev')`. For Linux, use `remotes::install_github("rspatial/terra")`. This package is also backwards compatible with the `raster` package. To install the development version of `raster` use `install.packages('raster', repos='https://rspatial.r-universe.dev')` on Windows or Mac or `remotes::install_github("rspatial/raster")` on Linux.
 
 ## Main Functions
 

README.md

@@ -1,12 +1,10 @@
 README
 ================
 Alexander Ilich
-January 30, 2022
+February 08, 2022
 
 # MultiscaleDEM
 
-## THIS IS A PRE-RELEASE. IT IS STILL BEING TESTED AND FUNCTIONS MAY RECEIVE MAJOR CHANGES.
-
 [![DOI](https://zenodo.org/badge/353158828.svg)](https://zenodo.org/badge/latestdoi/353158828)
 
 Please cite as
@@ -26,16 +24,21 @@ If you don’t already have remotes installed, use the code
 `install.packages("remotes")`
 
 Then to install this package use the code
-`remotes::install_github("ailich/MultiscaleDEM")` (you may need to
-install Rtools using the instructions found
+`remotes::install_github("ailich/MultiscaleDEM")` (If you are using
+Windows, you may need to install Rtools using the instructions found
 [here](https://cran.r-project.org/bin/windows/Rtools/)).
 
 This package relies on the `terra` package for handling of spatial
-raster data. To install the development version of `terra`, use
+raster data. The CRAN version of `terra` may not have all necessary
+functions, so it is recommended to install the development version of
+`terra`. The easiest way to install the development version of `terra`
+on Windows or Mac is to use
 `install.packages('terra', repos='https://rspatial.r-universe.dev')`.
-This package is also backwards compatible with the `raster` package. To
-install the development version of `raster` use
-`install.packages('raster', repos='https://rspatial.r-universe.dev')`
+For Linux, use `remotes::install_github("rspatial/terra")`. This package
+is also backwards compatible with the `raster` package. To install the
+development version of `raster` use
+`install.packages('raster', repos='https://rspatial.r-universe.dev')` on
+Windows or Mac or `remotes::install_github("rspatial/raster")` on Linux.
 
 ## Main Functions
 

inst/Terrain_Attributes_Explorer_App/app.R

@@ -46,13 +46,13 @@ classify_features2<- function(slope, planc, maxc, minc) {
     dplyr::case_when(is.na(slope) ~ NA_character_,
                  (slope > slope_tolerance) & (planc > curvature_tolerance) ~ "Ridge",
                  (slope > slope_tolerance) & (planc < -curvature_tolerance) ~ "Channel",
-                 slope > slope_tolerance ~ "Planar", #Maybe make this 7 (slope, new category)
+                 slope > slope_tolerance ~ "Planar Slope",
                  (maxc > curvature_tolerance) & (minc > curvature_tolerance) ~ "Peak",
                  (maxc > curvature_tolerance) & (minc < -curvature_tolerance) ~ "Pass",
                  maxc > curvature_tolerance ~ "Ridge",
                  (minc < -curvature_tolerance) & (maxc < -curvature_tolerance) ~ "Pit",
                  minc < -curvature_tolerance ~ "Channel",
-                 TRUE ~ "Planar")
+                 TRUE ~ "Planar Flat")
     }
 
 nr<-3