diff --git a/DESCRIPTION b/DESCRIPTION index 03c549c..4a951c3 100644 --- a/DESCRIPTION +++ b/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", diff --git a/R/AdjSD.R b/R/AdjSD.R index cf896f0..8dc04aa 100644 --- a/R/AdjSD.R +++ b/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)){ diff --git a/R/RDMV.R b/R/RDMV.R index 278b87a..a0f97a9 100644 --- a/R/RDMV.R +++ b/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 diff --git a/R/SAPA.R b/R/SAPA.R index 84e1a3f..9386e57 100644 --- a/R/SAPA.R +++ b/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)){ diff --git a/R/SlpAsp.R b/R/SlpAsp.R index 52d2165..287c3c9 100644 --- a/R/SlpAsp.R +++ b/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,27 +171,21 @@ 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) } @@ -199,15 +193,15 @@ SlpAsp <- function(r, w=c(3,3), unit="degrees", method="queen", metrics= c("slop 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) } diff --git a/R/SurfaceArea.R b/R/SurfaceArea.R index 617ee52..e0eef19 100644 --- a/R/SurfaceArea.R +++ b/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)){ diff --git a/R/TPI.R b/R/TPI.R index e2d1f30..0286d1c 100644 --- a/R/TPI.R +++ b/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 diff --git a/R/VRM.R b/R/VRM.R index 6408555..9c8684f 100644 --- a/R/VRM.R +++ b/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 diff --git a/R/WoodEvans.R b/R/WoodEvans.R index 5588e80..d8f80f5 100644 --- a/R/WoodEvans.R +++ b/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) } diff --git a/R/curvatures.R b/R/curvatures.R index ceb4b8a..62e6615 100644 --- a/R/curvatures.R +++ b/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) diff --git a/README.Rmd b/README.Rmd index aaa7570..e438fc2 100644 --- a/README.Rmd +++ b/README.Rmd @@ -37,13 +37,13 @@ This package relies on the `terra` package for handling of spatial raster data. ### Slope, Aspect and Curvature -- `SlpAsp` 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). +- `SlpAsp` calculates multi-scale slope and aspect according to 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 (Evans, 1980; Wood, 1996; Wilson et al. 2007). This is similar to [r.param.scale in GRASS GIS](https://grass.osgeo.org/grass80/manuals/r.param.scale.html). +- `WoodEvans` calculates slope, aspect, curvature, and morphometric features by fitting a quadratic surface to the focal window using ordinary least squares (Evans, 1980; Wilson et al., 2007; Wood, 1996). The morphometric features algorithm has been modified to use more robust measures of curvature based on the suggestions of Minár et al. (2020). ### Rugosity -- `VRM` - Vector ruggedness measure (Sappington et al. 2007). +- `VRM` - Vector ruggedness measure (Sappington et al. 2007) quantifies terrain ruggedness by measuring the dispersion of vectors orthogonal to the terrain surface. - `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. @@ -55,15 +55,15 @@ This package relies on the `terra` package for handling of spatial raster data. ### Relative Position -- `TPI` - Topographic Position Index (Weiss, 2001) +- `TPI` - Topographic Position Index (Weiss, 2001) is the difference between the value of a focal cell and the mean of the surrounding cells. -- `RDMV` - Relative Difference from Mean Value (Lecours et al., 2017) +- `RDMV` - Relative Difference from Mean Value (Lecours et al., 2017) is the difference between the value of a focal cell and the mean of the cells in the focal window divided by the range of the values in the focal window. -- `BPI` - Bathymetric Position Index (Lundblad et al., 2006) +- `BPI` - Bathymetric Position Index (Lundblad et al., 2006) is the difference between the value of a focal cell and the mean of the surrounding cells contained within an annulus shaped window. ## 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. The only metric that does not follow this syntax is BPI which uses an annulus shaped focal window. +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. The only metric that does not follow this syntax is BPI which uses an annulus shaped focal window which we will calculate using an inner radius of 2 and an outer radius of 4 cells. **Load packages** @@ -225,11 +225,6 @@ tm_shape(rdmv, raster.downsample = FALSE)+ tm_layout(legend.outside = TRUE, main.title="RDMV") ``` -BPI is a modification of TPI that uses an annulus shaped focal window and therefore requires an inner and outer radius. This can be specified in cell units (number of raster cells) or in map units (e.g. meters) which can be useful if your x and y resolutions are not equal. For example, an annulus window with an inner radius of 2 cells and an outer radius of 4 cells would be -```{r} -annulus_window(radius = c(2,4), unit = "cell") -``` - ```{r} bpi<- BPI(r, radius = c(2,4), unit = "cell", na.rm = TRUE) ``` @@ -240,6 +235,12 @@ tm_shape(bpi, raster.downsample = FALSE)+ tm_layout(legend.outside = TRUE, main.title="BPI") ``` +BPI is a modification of TPI that uses an annulus shaped focal window and therefore requires an inner and outer radius. This can be specified in cell units (number of raster cells) or in map units (e.g. meters) which can be useful if your x and y resolutions are not equal. For example, an annulus window with an inner radius of 2 cells and an outer radius of 4 cells would be + +```{r} +annulus_window(radius = c(2,4), unit = "cell") +``` + # 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 @@ -256,6 +257,8 @@ Lecours, V., Devillers, R., Simms, A.E., Lucieer, V.L., Brown, C.J., 2017. Towar Lundblad, E.R., Wright, D.J., Miller, J., Larkin, E.M., Rinehart, R., Naar, D.F., Donahue, B.T., Anderson, S.M., Battista, T., 2006. A benthic terrain classification scheme for American Samoa. Marine Geodesy 29, 89–111. +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 + 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. diff --git a/README.md b/README.md index 8986d8d..c7631c6 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ README ================ Alexander Ilich -December 08, 2021 +December 09, 2021 # MultiscaleDEM @@ -41,20 +41,23 @@ install the development version of `raster` use ### Slope, Aspect and Curvature -- `SlpAsp` 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). +- `SlpAsp` calculates multi-scale slope and aspect according to 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 (Evans, 1980; Wood, 1996; Wilson et - al. 2007). This is similar to [r.param.scale in GRASS - GIS](https://grass.osgeo.org/grass80/manuals/r.param.scale.html). + ordinary least squares (Evans, 1980; Wilson et al., 2007; Wood, + 1996). The morphometric features algorithm has been modified to use + more robust measures of curvature based on the suggestions of Minár + et al. (2020). ### Rugosity -- `VRM` - Vector ruggedness measure (Sappington et al. 2007). +- `VRM` - Vector ruggedness measure (Sappington et al. 2007) + quantifies terrain ruggedness by measuring the dispersion of vectors + orthogonal to the terrain surface. - `SAPA` - Calculates the Surface Area to Planar Area (Jenness, 2004). Additionally, planar area can be corrected for slope (Du Preez @@ -76,11 +79,18 @@ install the development version of `raster` use ### Relative Position -- `TPI` - Topographic Position Index (Weiss, 2001) +- `TPI` - Topographic Position Index (Weiss, 2001) is the difference + between the value of a focal cell and the mean of the surrounding + cells. - `RDMV` - Relative Difference from Mean Value (Lecours et al., 2017) + is the difference between the value of a focal cell and the mean of + the cells in the focal window divided by the range of the values in + the focal window. -- `BPI` - Bathymetric Position Index (Lundblad et al., 2006) +- `BPI` - Bathymetric Position Index (Lundblad et al., 2006) is the + difference between the value of a focal cell and the mean of the + surrounding cells contained within an annulus shaped window. ## Tutorial @@ -90,7 +100,8 @@ 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. The only metric that does not follow this syntax is BPI which uses an annulus shaped focal -window. +window which we will calculate using an inner radius of 2 and an outer +radius of 4 cells. **Load packages** @@ -165,6 +176,12 @@ rdmv<- RDMV(r, w=c(5,5), na.rm = TRUE) ![](README_files/figure-gfm/unnamed-chunk-18-1.png) +``` r +bpi<- BPI(r, radius = c(2,4), unit = "cell", na.rm = TRUE) +``` + +![](README_files/figure-gfm/unnamed-chunk-20-1.png) + BPI is a modification of TPI that uses an annulus shaped focal window and therefore requires an inner and outer radius. This can be specified in cell units (number of raster cells) or in map units (e.g. meters) @@ -187,12 +204,6 @@ annulus_window(radius = c(2,4), unit = "cell") ## [8,] NA NA 1 1 1 1 1 NA NA ## [9,] NA NA NA NA 1 NA NA NA NA -``` r -bpi<- BPI(r, radius = c(2,4), unit = "cell", na.rm = TRUE) -``` - -![](README_files/figure-gfm/unnamed-chunk-21-1.png) - # References Du Preez, C., 2015. A new arc–chord ratio (ACR) rugosity index for @@ -224,6 +235,12 @@ Naar, D.F., Donahue, B.T., Anderson, S.M., Battista, T., 2006. A benthic terrain classification scheme for American Samoa. Marine Geodesy 29, 89–111. +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. + + 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. diff --git a/README_files/figure-gfm/unnamed-chunk-20-1.png b/README_files/figure-gfm/unnamed-chunk-20-1.png new file mode 100644 index 0000000..a5057bb Binary files /dev/null and b/README_files/figure-gfm/unnamed-chunk-20-1.png differ diff --git a/man/AdjSD.Rd b/man/AdjSD.Rd index 69dea53..53fbdfc 100644 --- a/man/AdjSD.Rd +++ b/man/AdjSD.Rd @@ -11,7 +11,7 @@ AdjSD(r, w = c(3, 3), na.rm = FALSE, include_scale = FALSE) \item{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.} -\item{na.rm}{A logical vector indicating whether or not to remove NA values before calculations} +\item{na.rm}{A logical indicating whether or not to remove NA values before calculations} \item{include_scale}{logical indicating whether to append window size to the layer names (default = FALSE)} } diff --git a/man/RDMV.Rd b/man/RDMV.Rd index d40c882..c5ddf45 100644 --- a/man/RDMV.Rd +++ b/man/RDMV.Rd @@ -11,7 +11,7 @@ RDMV(r, w = c(3, 3), na.rm = FALSE, pad = FALSE, include_scale = FALSE) \item{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.} -\item{na.rm}{A logical vector indicating whether or not to remove NA values before calculations} +\item{na.rm}{A logical indicating whether or not to remove NA values before calculations} \item{include_scale}{logical indicating whether to append window size to the layer names (default = FALSE)} } diff --git a/man/TPI.Rd b/man/TPI.Rd index 55923bd..2ff3ae7 100644 --- a/man/TPI.Rd +++ b/man/TPI.Rd @@ -11,7 +11,7 @@ TPI(r, w = c(3, 3), na.rm = FALSE, include_scale = FALSE) \item{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.} -\item{na.rm}{A logical vector indicating whether or not to remove NA values before calculations} +\item{na.rm}{A logical indicating whether or not to remove NA values before calculations} \item{include_scale}{logical indicating whether to append window size to the layer names (default = FALSE)} } diff --git a/man/VRM.Rd b/man/VRM.Rd index 23832e1..906c59d 100644 --- a/man/VRM.Rd +++ b/man/VRM.Rd @@ -11,7 +11,7 @@ VRM(r, w, na.rm = FALSE, include_scale = FALSE) \item{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.} -\item{na.rm}{A logical vector indicating whether or not to remove NA values before calculations} +\item{na.rm}{A logical indicating whether or not to remove NA values before calculations} \item{include_scale}{logical indicating whether to append window size to the layer names (default = FALSE)} } diff --git a/man/WoodEvans.Rd b/man/WoodEvans.Rd index 86bc46a..c9707a9 100644 --- a/man/WoodEvans.Rd +++ b/man/WoodEvans.Rd @@ -31,7 +31,7 @@ WoodEvans( \item{curvature_tolerance}{Curvature tolerance that defines 'planar' surface (default = 0.0001). Relevant for the features layer.} -\item{na.rm}{Logical vector indicating whether or not to remove NA values before calculations.} +\item{na.rm}{Logical indicating whether or not to remove NA values before calculations.} \item{include_scale}{Logical indicating whether to append window size to the layer names (default = FALSE).} diff --git a/man/kmax.Rd b/man/kmax.Rd new file mode 100644 index 0000000..2e11ff1 --- /dev/null +++ b/man/kmax.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{kmax} +\alias{kmax} +\title{Calculate max curvature} +\usage{ +kmax(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +Calculate max curvature, kmax, from the equation Z =ax^2+by^2+cxy+dx+e+f. +} +\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 +} diff --git a/man/kmean.Rd b/man/kmean.Rd new file mode 100644 index 0000000..6c8772c --- /dev/null +++ b/man/kmean.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{kmean} +\alias{kmean} +\title{Calculate mean curvature} +\usage{ +kmean(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +Calculate mean curvature, kmean, from the equation Z =ax^2+by^2+cxy+dx+e+f. +} +\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 +} diff --git a/man/kmin.Rd b/man/kmin.Rd new file mode 100644 index 0000000..f8987dd --- /dev/null +++ b/man/kmin.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{kmin} +\alias{kmin} +\title{Calculate min curvature} +\usage{ +kmin(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +Calculate min curvature, kmin, from the equation Z =ax^2+by^2+cxy+dx+e+f. +} +\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 +} diff --git a/man/knc.Rd b/man/knc.Rd new file mode 100644 index 0000000..2966ea6 --- /dev/null +++ b/man/knc.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{knc} +\alias{knc} +\title{Calculate normal contour curvature} +\usage{ +knc(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +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. +} +\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 +} diff --git a/man/kns.Rd b/man/kns.Rd new file mode 100644 index 0000000..e22fe66 --- /dev/null +++ b/man/kns.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{kns} +\alias{kns} +\title{Calculate normal slope line curvature} +\usage{ +kns(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +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. +} +\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 +} diff --git a/man/ku.Rd b/man/ku.Rd new file mode 100644 index 0000000..0f4aad1 --- /dev/null +++ b/man/ku.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{ku} +\alias{ku} +\title{Calculate unsphericity curvature} +\usage{ +ku(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +Calculate unsphericity curvature, ku, from the equation Z =ax^2+by^2+cxy+dx+e+f. +} +\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 +} diff --git a/man/tgc.Rd b/man/tgc.Rd new file mode 100644 index 0000000..f038fc6 --- /dev/null +++ b/man/tgc.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/curvatures.R +\name{tgc} +\alias{tgc} +\title{Calculate contour geodesic torsion} +\usage{ +tgc(a, b, c, d, e, f) +} +\arguments{ +\item{a}{regression coefficient} + +\item{b}{regression coefficient} + +\item{c}{regression coefficient} + +\item{d}{regression coefficient} + +\item{e}{regression coefficient} + +\item{f}{regression intercept} +} +\description{ +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. +} +\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 +}