Updates for pre-release

- Add TPI and RDMV
- Add ReadME
- Add quotes in messages for SlopeAspect
- Add authors in Description

DESCRIPTION

@@ -1,15 +1,24 @@
 Package: MultiscaleDEM
-Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9000
+Title: Calculates multi-scale geomorphometric terrain attributes from regularly gridded DEM/bathymetry rasters
+Version: 0.0.1
 Authors@R: 
-    person(given = "First",
-           family = "Last",
+     c(
+    person(given = "Alexander",
+           family = "Ilich",
            role = c("aut", "cre"),
-           email = "[email protected]",
-           comment = c(ORCID = "YOUR-ORCID-ID"))
-Description: What the package does (one paragraph).
-License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
-    license
+           email = "[email protected]",
+           comment = c(ORCID = "0000-0003-1758-8499")),
+    person(given = "Vincent",
+           family = "Lecours",
+           role = c("aut")),
+    person(given = "Benjamin",
+           family = "Misiuk",
+           role = c("aut")),
+    person(given = "Steven",
+           family = "Murawski",
+           role = c("aut")))
+Description: Calculates multi-scale geomorphometric terrain attributes from regularly gridded DEM/bathymetry rasters
+License: `use_mit_license()`
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -1,9 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 export(AdjSD)
+export(RDMV)
 export(SAPA)
 export(SlopeAspect)
 export(SurfaceArea)
+export(TPI)
 export(VRM)
 export(WoodEvans)
 import(raster)

R/RDMV.R

@@ -0,0 +1,28 @@
+#' Calculates Relative Difference from Mean Value (RDMV)
+#'
+#' Calculates Relative Difference from Mean Value (RDMV). RDMV = (focal_value - local_mean)/local_range
+#' @param r DEM as a raster layer
+#' @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 pad logical value specifying whether rows/columns of NA's should be padded to the edge of the raster to remove edge effects (FALSE by default). If pad is TRUE, na.rm must be TRUE.
+#' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
+#' @return a RasterLayer
+#' @import raster
+#' @export
+
+RDMV<- function(r, w=c(3,3), na.rm=FALSE, pad=FALSE, include_scale=FALSE){
+  if(length(w)==1){w<- rep(w, 2)}
+  if(any(w<3) | any(0 == (w %% 2))){
+    stop("Error: w must be odd and greater than or equal to 3")}
+
+  w_mat <- matrix(1, nrow=w[1], ncol=w[2])
+
+  localmean<- raster::focal(x = r, w= w_mat, fun=mean, na.rm = na.rm, pad=pad)
+  localmax<- raster::focal(x = r, w= w_mat, fun=max, na.rm = na.rm, pad=pad)
+  localmin<- raster::focal(x = r, w= w_mat, fun=min, na.rm = na.rm, pad=pad)
+  rdmv<- (r - localmean)/(localmax-localmin)
+  names(rdmv)<- "rdmv"
+
+  if(include_scale){names(rdmv)<- paste0(names(rdmv), "_", w[1],"x", w[2])} #Add scale to layer names
+  return(rdmv)
+}

R/SlopeAspect.R

@@ -8,6 +8,7 @@
 #' @param metrics a character string or vector of character strings of which terrain atrributes to return ("slope" and/or "aspect"). Default is c("slope", "aspect").
 #' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
 #' @param mask_aspect A logical. If slope evaluates to 0, aspect will be set to NA when mask_aspect is TRUE (the default). If FALSE, when slope is 0 aspect will be pi/2 radians or 90 degrees which is the behavior of raster::terrain.
+#' @param mask_aspect A logical. If slope evaluates to 0, aspect will be set to NA when mask_aspect is TRUE (the default). If FALSE, when slope is 0 aspect will be pi/2 radians or 90 degrees which is the behavior of raster::terrain.
 #' @return a RasterStack or RasterLayer of slope and/or aspect
 #' @details When method="rook", slope and aspect are computed according to Fleming and Hoffer (1979) and Ritter (1987). When method="queen", slope and aspect are computed according to Horn (1981). These are the standard slope algorithms found in many GIS packages but are traditionally restricted to a 3 x 3 window size. Misiuk et al (2021) extended these classical formulations  to multiple window sizes. This function modifies the code from Misiuk et al (2021) to allow for rectangular rather than only square windows and also added aspect.
 #' @import raster
@@ -22,10 +23,10 @@ SlopeAspect <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c(
     stop("w must be odd")
   }
   if(!(unit %in% c("degrees", "radians"))){
-    stop("unit must be degrees or radians")
+    stop("unit must be `degrees` or `radians`")
   }
   if(!(method %in% c("queen", "rook"))){
-    stop("method must be queen or rook")
+    stop("method must be `queen` or `rook`")
   }
 
   if(any(!(metrics %in% c("slope", "aspect","northness", "eastness")))){

R/TPI.R

@@ -0,0 +1,35 @@
+#' Calculates Topographic Position Index
+#'
+#' 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 raster layer
+#' @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 pad logical value specifying whether rows/columns of NA's should be padded to the edge of the raster to remove edge effects (FALSE by default). If pad is TRUE, na.rm must be TRUE.
+#' @param include_scale logical indicating whether to append window size to the layer names (default = FALSE)
+#' @return a RasterLayer
+#' @import raster
+#' @export
+#' 
+
+TPI<- function(r, w=c(3,3), na.rm=FALSE, pad=FALSE, include_scale=FALSE){
+  if(length(w)==1){w<- rep(w, 2)}
+
+  if(any(w<3) | any(0 == (w %% 2))){
+    stop("Error: w must be odd and greater than or equal to 3")}
+
+  w_mat <- matrix(1, nrow=w[1], ncol=w[2])
+  center_idx<- ceiling(0.5 * length(w_mat))
+  w_mat[center_idx] <- 0 
+  if(!na.rm){
+    w_mat[w_mat>0]<- 1/sum(w_mat>0)
+    tpi<- r - raster::focal(x=r, w=w_mat, fun=sum, na.rm=FALSE, pad=FALSE)
+    } else{
+      custom_mean_fun<- function(x, na.rm) {
+        return(mean(x[-center_idx], na.rm=na.rm))
+        } #Using terra::focal with NA weights matrix will be better
+      tpi<- r - raster::focal(x = r, w = w_mat, fun = custom_mean_fun, na.rm = TRUE, pad = pad)
+    }
+  names(tpi)<- "TPI"
+  if(include_scale){names(tpi)<- paste0(names(tpi), "_", w[1],"x", w[2])} #Add scale to layer names
+return(tpi)
+}

README.Rmd

@@ -0,0 +1,237 @@
+---
+title: "README"
+author: "Alexander Ilich"
+date: "`r format(Sys.time(), '%B %d, %Y')`"
+output: 
+  github_document:
+    pandoc_args: --webtex
+---
+# MultiscaleDEM
+
+## THIS IS A PRE-RELEASE. IT IS STILL BEING TESTED AND FUNCTIONS MAY RECEIVE MAJOR CHANGES.
+
+DOI GOES HERE
+
+Please cite as
+
+CITATION GOES HERE
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) #Default chunk options
+```
+
+## Purpose
+
+This package calculates multi-scale geomorphometric terrain attributes from regularly gridded DEM/bathymetry rasters.
+
+## Install and Load Package
+
+If you don't already have devtools installed, use the code `install.packages("devtools")`
+
+Then to install this package use the code `devtools::install_github("ailich/MultiscaleDEM")`
+
+## Main Functions
+
+### Slope, Aspect and Curvature
+
+-   `SlopeAspect` calculates multi-scale slope and aspect according to the Misiuk et al (2021) which is a modification of the traditional 3 x 3 slope and aspect algorithms (Fleming and Hoffer, 1979; Horn et al., 1981; Ritter, 1987).
+
+-   `WoodEvans` calculates slope, aspect, curvature, and morphometric features by fitting a quadratic surface to the focal window using ordinary least squares (Wood, 1996). This is similar to [r.param.scale in GRASS GIS](https://grass.osgeo.org/grass80/manuals/r.param.scale.html).
+
+### Rugosity
+
+-   `VRM` - Vector ruggedness measure (Sappington et al. 2007).
+
+-   `SAPA` - Calculates the Surface Area to Planar Area (Jenness, 2004). Additionally, planar area can be corrected for slope (Du Preez 2015). Additionally, a proposed extension to multiple scales is provided by summing the surface areas within the focal window and adjusting the planar area of the focal window using multi-scale slope.
+
+-   `SurfaceArea` - Calculate the surface area of each grid cell (Jenness, 2004). Used within `SAPA`.
+
+-   `AdjSD`- This new proposed rugosity metric modifies the standard deviation of elevation/bathymetry to account for slope. It does this by first fitting a plane to the data in the focal window using ordinary least squares, and then extracting the residuals, and then calculating the standard deviation of the residuals.
+
+![](images/adj_sd.png)
+
+### Relative Position
+
+-   `TPI` - Topographic Position Index (Weiss, 2001)
+
+-   `RDMV` - Relative Difference from Mean Value (Lecours et al., 2017)
+
+## Tutorial
+
+In this tutorial we will calculate various terrain attributes using a 5 x 5 cell rectangular window. Any rectangular odd numbered window size however could be used. Window sizes are specified with a vector of length 2 of `c(n_rows, n_cols)`. If a single number is provided it will be used for both the number of rows and columns.
+
+Load packages
+
+```{r message=FALSE}
+library(raster) #Load raster package
+library(MultiscaleDEM) #Load MultiscaleDEM package
+```
+
+See package help page
+
+```{r eval=FALSE}
+help(package="MultiscaleDEM")
+```
+
+Read in Data
+```{r}
+r<- raster(volcano, xmn=0, xmx=ncol(volcano)*10, ymn=0, ymx=nrow(volcano)*10) #Use built in volcano data set
+crs(r)<- CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs") #Give r a projection. This is not the right projection, but some functions currently throw an error if there is no crs.
+```
+
+```{r echo= FALSE, message=FALSE}
+library(tmap) #For plotting
+tm_shape(r, raster.downsample = FALSE)+
+  tm_raster(palette = colorRamps::matlab.like(100), style = "cont", legend.reverse = TRUE, title = "")+
+  tm_layout(legend.outside=TRUE, main.title= "Elevation")
+```
+
+### Slope, Aspect, and Curvature
+
+```{r}
+slp_asp<- SlopeAspect(r = r, w = c(5,5), unit = "degrees", method = "queen")
+```
+
+```{r echo=FALSE}
+slp_asp_list<- vector(mode="list", length = nlayers(slp_asp))
+for (i in 1:length(slp_asp_list)) {
+  curr_var<- names(slp_asp)[i]
+  
+  if (grepl(pattern = "(^northness)|(eastness)", curr_var)) {
+    breaks<- c(-1,0,1)
+    midpoint<- 0
+    curr_pal<- c("blue", "gray", "red")
+  } else if (grepl(pattern = "aspect", curr_var)) {
+    curr_pal<-c("blue", "purple", "red", "orange", "yellow", "green", "cyan", "blue")
+    breaks<- c(0,90,180,270,360)
+    midpoint<- 180
+  } else{
+    curr_pal<- colorRamps::matlab.like(100)
+    midpoint<- NULL
+    breaks<- NULL
+  }
+  
+  slp_asp_list[[i]]<- tm_shape(slp_asp[[i]], raster.downsample = FALSE) +
+    tm_raster(palette = curr_pal, style= "cont", title = "", breaks = breaks, midpoint = midpoint, legend.reverse = TRUE)+
+      tm_layout(main.title = curr_var, 
+      main.title.position = "center",
+      main.title.size=0.75)
+  }
+slp_asp_plot<- tmap_arrange(slp_asp_list, ncol=2)
+slp_asp_plot
+```
+
+```{r}
+WE<- WoodEvans(r, w = c(5,5), unit = "degrees", return_aspect = TRUE, na.rm = TRUE, pad = TRUE)
+```
+
+```{r echo= FALSE}
+WE_list<- vector(mode="list", length = nlayers(WE))
+for (i in 1:length(WE_list)) {
+  curr_var<- names(WE)[i]
+  
+  if (grepl(pattern = "(^Northness)|(Eastness)", curr_var)) {
+    breaks<- c(-1,0,1)
+    midpoint<- 0
+    curr_pal<- c("blue", "gray", "red")
+    style<- "cont"
+  } else if (grepl(pattern = "Aspect", curr_var)) {
+    curr_pal<-c("blue", "purple", "red", "orange", "yellow", "green", "cyan", "blue")
+    breaks<- c(0,90,180,270,360)
+    midpoint<- 180
+    style<- "cont"
+    } else if(grepl(pattern = "^Features", curr_var)) {
+    curr_pal<- c("gray", "black", "blue", "green", "yellow", "red")
+    midpoint<- NULL
+    breaks<- NULL
+    style<- "cat"
+  } else if(grepl(pattern = "Curv", curr_var)){
+    curr_pal<- c("blue", "gray", "red")
+    style<- "cont"
+  } else{
+    curr_pal<- colorRamps::matlab.like(100)
+    midpoint<- NULL
+    breaks<- NULL
+    style<- "cont"}
+  
+  WE_list[[i]]<- tm_shape(WE[[i]], raster.downsample = FALSE) +
+    tm_raster(palette = curr_pal, style= style, title = "", breaks = breaks, midpoint = midpoint, legend.reverse = TRUE)+
+      tm_layout(main.title = curr_var, 
+      main.title.position = "center",
+      main.title.size=0.75)
+}
+WE_plot<- tmap_arrange(WE_list, ncol=4)
+WE_plot
+```
+### Rugosity
+
+```{r}
+vrm<- VRM(r, w=c(5,5))
+```
+
+```{r echo=FALSE}
+tm_shape(vrm, raster.downsample = FALSE)+
+  tm_raster(palette = colorRamps::matlab.like(100), style = "cont", legend.reverse = TRUE, title="")+
+  tm_layout(legend.outside = TRUE, main.title="VRM")
+```
+
+Note: multi-scale SAPA is experimental. The established metric by De Preez (2015) would use `w=1`.
+```{r}
+sapa<- SAPA(r, w=c(5,5), slope_correction = TRUE)
+```
+
+```{r echo=FALSE}
+tm_shape(sapa, raster.downsample = FALSE)+
+  tm_raster(palette = colorRamps::matlab.like(100), style = "cont", legend.reverse = TRUE, title="")+
+  tm_layout(legend.outside = TRUE, main.title="SAPA")
+```
+```{r}
+adj_SD<- AdjSD(r, w=c(5,5), na.rm = TRUE, pad=TRUE)
+```
+
+```{r echo=FALSE}
+tm_shape(adj_SD, raster.downsample = FALSE)+
+  tm_raster(palette = colorRamps::matlab.like(100), style = "cont", legend.reverse = TRUE, title="")+
+  tm_layout(legend.outside = TRUE, main.title="Adjusted SD")
+```
+### Relative Position
+```{r}
+tpi<- TPI(r, w=c(5,5), na.rm = TRUE, pad = TRUE)
+```
+
+```{r echo=FALSE}
+tm_shape(tpi, raster.downsample = FALSE)+
+  tm_raster(palette = c("blue", "gray", "red"), style = "cont", midpoint=0, legend.reverse = TRUE, title="")+
+  tm_layout(legend.outside = TRUE, main.title="TPI")
+```
+```{r}
+rdmv<- RDMV(r, w=c(5,5), na.rm = TRUE, pad = TRUE)
+```
+
+```{r echo=FALSE}
+tm_shape(rdmv, raster.downsample = FALSE)+
+  tm_raster(palette = c("blue", "gray", "red"), style = "cont", midpoint=0, legend.reverse = TRUE, title="")+
+  tm_layout(legend.outside = TRUE, main.title="RDMV")
+```
+
+# 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
+
+Fleming, M.D., Hoffer, R.M., 1979. Machine processing of landsat MSS data and DMA topographic data for forest cover type mapping (No. LARS Technical Report 062879). Laboratory for Applications of Remote Sensing, Purdue University, West Lafayette, Indiana.
+
+Horn, B.K., 1981. Hill Shading and the Reflectance Map. Proceedings of the IEEE 69, 14–47.
+
+Jenness, J.S., 2004. Calculating landscape surface area from digital elevation models. Wildlife Society Bulletin 32, 829–839. https://doi.org/10.2193/0091-7648(2004)032[0829:CLSAFD]2.0.CO;2
+
+Lecours, V., Devillers, R., Simms, A.E., Lucieer, V.L., Brown, C.J., 2017. Towards a Framework for Terrain Attribute Selection in Environmental Studies. Environmental Modelling & Software 89, 19–30. https://doi.org/10.1016/j.envsoft.2016.11.027
+
+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
+
+Ritter, P., 1987. A vector-based slope and aspect generation algorithm. Photogrammetric Engineering and Remote Sensing 53, 1109–1111.
+
+Weiss, A., 2001. Topographic Position and Landforms Analysis. Presented at the ESRI user conference, San Diego, CA.
+
+Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.
+
+