diff --git a/DESCRIPTION b/DESCRIPTION index 4b2babd..b368cdf 100644 --- a/DESCRIPTION +++ b/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", diff --git a/NAMESPACE b/NAMESPACE index fef0dc2..90ae9c1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(AdjSD) export(BPI) export(DMV) export(DirSlp) +export(Pfit) export(Qfit) export(RIE) export(RelPos) diff --git a/R/AdjSD.R b/R/AdjSD.R index 18994f2..2fd34d5 100644 --- a/R/AdjSD.R +++ b/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 diff --git a/R/BPI.R b/R/BPI.R index 2289f02..57d9540 100644 --- a/R/BPI.R +++ b/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 diff --git a/R/DMV.R b/R/DMV.R index d174553..e47e934 100644 --- a/R/DMV.R +++ b/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 diff --git a/R/DirSlope.R b/R/DirSlope.R index 5edf131..af86ae8 100644 --- a/R/DirSlope.R +++ b/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 diff --git a/R/Pfit.R b/R/Pfit.R index 6e57b84..6409275 100644 --- a/R/Pfit.R +++ b/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,21 +63,51 @@ 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) @@ -73,12 +115,89 @@ Pfit<- function(r, w=c(3,3), unit= "degrees",na.rm=FALSE, include_scale=FALSE, f } # 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) } diff --git a/R/Qfit.R b/R/Qfit.R index 6304384..8c1cbee 100644 --- a/R/Qfit.R +++ b/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){ diff --git a/R/RIE.R b/R/RIE.R index 1c72c78..f002ab2 100644 --- a/R/RIE.R +++ b/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 diff --git a/R/RcppExports.R b/R/RcppExports.R index d4ac439..5075e9a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,32 +1,24 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -subset_mat_rows <- function(Input_Matrix, Input_Log_Vec) { - .Call(`_MultiscaleDTM_subset_mat_rows`, Input_Matrix, Input_Log_Vec) -} - -C_OLS_params <- function(X, Y) { - .Call(`_MultiscaleDTM_C_OLS_params`, X, Y) -} - -C_OLS_params2 <- function(Xt, XtX_inv, Y) { - .Call(`_MultiscaleDTM_C_OLS_params2`, Xt, XtX_inv, Y) +C_AdjSD_narmT <- function(z, X_full, ni, nw) { + .Call(`_MultiscaleDTM_C_AdjSD_narmT`, z, X_full, ni, nw) } -C_OLS_resid <- function(X, Y) { - .Call(`_MultiscaleDTM_C_OLS_resid`, X, Y) +C_AdjSD_narmF <- function(z, X, Xt, XtX_inv, ni, nw) { + .Call(`_MultiscaleDTM_C_AdjSD_narmF`, z, X, Xt, XtX_inv, ni, nw) } -C_OLS_resid2 <- function(X, Xt, XtX_inv, Y) { - .Call(`_MultiscaleDTM_C_OLS_resid2`, X, Xt, XtX_inv, Y) +C_CountVals <- function(z, ni, nw) { + .Call(`_MultiscaleDTM_C_CountVals`, z, ni, nw) } -C_Qfit1_narmT <- function(z, X_full, ni, nw) { - .Call(`_MultiscaleDTM_C_Qfit1_narmT`, z, X_full, ni, nw) +C_Qfit1_narmT <- function(z, X_full, return_intercept, ni, nw) { + .Call(`_MultiscaleDTM_C_Qfit1_narmT`, z, X_full, return_intercept, ni, nw) } -C_Qfit1_narmF <- function(z, X, Xt, XtX_inv, ni, nw) { - .Call(`_MultiscaleDTM_C_Qfit1_narmF`, z, X, Xt, XtX_inv, ni, nw) +C_Qfit1_narmF <- function(z, X, Xt, XtX_inv, return_intercept, ni, nw) { + .Call(`_MultiscaleDTM_C_Qfit1_narmF`, z, X, Xt, XtX_inv, return_intercept, ni, nw) } C_Qfit2_narmT <- function(z, X_full, ni, nw) { @@ -37,18 +29,6 @@ C_Qfit2_narmF <- function(z, X, Xt, XtX_inv, ni, nw) { .Call(`_MultiscaleDTM_C_Qfit2_narmF`, z, X, Xt, XtX_inv, ni, nw) } -C_AdjSD_narmT <- function(z, X_full, ni, nw) { - .Call(`_MultiscaleDTM_C_AdjSD_narmT`, z, X_full, ni, nw) -} - -C_AdjSD_narmF <- function(z, X, Xt, XtX_inv, ni, nw) { - .Call(`_MultiscaleDTM_C_AdjSD_narmF`, z, X, Xt, XtX_inv, ni, nw) -} - -C_Pfit1_narmF <- function(z, X, Xt, XtX_inv, idx, ni, nw) { - .Call(`_MultiscaleDTM_C_Pfit1_narmF`, z, X, Xt, XtX_inv, idx, ni, nw) -} - C_TriArea <- function(a, b, c) { .Call(`_MultiscaleDTM_C_TriArea`, a, b, c) } @@ -57,7 +37,3 @@ C_SurfaceArea <- function(z, x_res, y_res, na_rm, ni, nw) { .Call(`_MultiscaleDTM_C_SurfaceArea`, z, x_res, y_res, na_rm, ni, nw) } -C_CountVals <- function(z, ni, nw) { - .Call(`_MultiscaleDTM_C_CountVals`, z, ni, nw) -} - diff --git a/R/RelPos.R b/R/RelPos.R index d114964..19315ba 100644 --- a/R/RelPos.R +++ b/R/RelPos.R @@ -110,9 +110,7 @@ annulus_window<- function(radius, unit, resolution){ #' @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() #' rpos<- RelPos(r, w = c(5,5), shape= "rectangle", exclude_center = TRUE, na.rm = TRUE) #' plot(rpos) #' @import terra diff --git a/R/SAPA.R b/R/SAPA.R index 8c4833f..cb1bb0f 100644 --- a/R/SAPA.R +++ b/R/SAPA.R @@ -14,9 +14,7 @@ #' @param overwrite logical. If TRUE, filename is overwritten (default is FALSE). #' @param wopt list with named options for writing files as in writeRaster #' @examples -#' r<- rast(volcano, extent= ext(2667400, 2667400 + -#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -#' crs = "EPSG:27200") +#' r<- erupt() #' sapa<- SAPA(r, w=c(5,5), slope_correction = TRUE) #' plot(sapa) #' @details Planar area is calculated as the x_dis * y_dis if uncorrected for slope and (x_dis * y_dis)/cos(slope) if corrected for slope. When w=1, this is called "native" scale and is equivalent to what is presented in Du Preez (2015) and available in the ArcGIS Benthic Terrain Modeller add-on. In this case operations are performed on a per cell basis where x_dis is the resolution of the raster in the x direction (left/right) and y_dis is the resolution of the raster in the y direction (up/down) and slope is calculated using the Horn (1981) method. To expand this to multiple scales of analysis, at w > 1 slope is calculated based on Misiuk et al (2021) which provides a modification of the Horn method to extend the matric to multiple spatial scales. Planar area is calculated the same way as for w=1 except that now x_dis is the x resolution of the raster * the number of columns in the focal window, and y_dis is y resolution of the raster * the number of rows. For w > 1, surface area is calculated as the sum of surface areas within the focal window. Although the (modified) Horn slope is used by default to be consistent with Du Preez (2015), slope calculated using a different algorithm (e.g. Wood 1996) could be supplied using the slope_layer argument. Additionally, a slope raster can be supplied if you have already calculated it and do not wish to recalculate it. However, be careful to supply a slope layer measured in radians and calculated at the relevant scale (2 larger than the w of SAPA). diff --git a/R/SlpAsp.R b/R/SlpAsp.R index 7c58fe7..02dd75d 100644 --- a/R/SlpAsp.R +++ b/R/SlpAsp.R @@ -18,9 +18,7 @@ #' @importFrom raster writeRaster #' @return a SpatRaster or RasterStack of slope and/or aspect (and components of aspect) #' @examples -#' r<- rast(volcano, extent= ext(2667400, 2667400 + -#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -#' crs = "EPSG:27200") +#' r<- erupt() #' slp_asp<- SlpAsp(r = r, w = c(5,5), unit = "degrees", #' method = "queen", metrics = c("slope", "aspect", #' "eastness", "northness")) diff --git a/R/SurfaceArea.R b/R/SurfaceArea.R index 2b62c80..8050974 100644 --- a/R/SurfaceArea.R +++ b/R/SurfaceArea.R @@ -8,9 +8,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() #' sa<- SurfaceArea(r) #' plot(sa) #' @import terra diff --git a/R/TPI.R b/R/TPI.R index c25c804..897303f 100644 --- a/R/TPI.R +++ b/R/TPI.R @@ -13,9 +13,7 @@ #' @param wopt List with named options for writing files as in writeRaster. #' @return SpatRaster or RasterLayer. #' @examples -#' r<- rast(volcano, extent= ext(2667400, 2667400 + -#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -#' crs = "EPSG:27200") +#' r<- erupt() #' tpi<- TPI(r, w=c(5,5), shape="rectangle", stand="none", na.rm = TRUE) #' plot(tpi) #' @import terra diff --git a/R/VRM.R b/R/VRM.R index 46a4ab9..e9131e2 100644 --- a/R/VRM.R +++ b/R/VRM.R @@ -10,9 +10,7 @@ #' @param wopt list with named options for writing files as in writeRaster #' @return a RasterLayer #' @examples -#' r<- rast(volcano, extent= ext(2667400, 2667400 + -#' ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -#' crs = "EPSG:27200") +#' r<- erupt() #' vrm<- VRM(r, w=c(5,5), na.rm = TRUE) #' plot(vrm) #' @import terra @@ -62,23 +60,37 @@ VRM<- function(r, w=c(3,3), na.rm = FALSE, include_scale=FALSE, filename=NULL, o } #SlpAsp can use na.rm in slope calculations, terra::terrain cannot but can handle spherical geometry. #Decompose vectors into components - sin.slp <- terra::math(sa$slope, fun="sin", wopt=wopt) - Xrast <- terra::math(sa$aspect, fun="sin", wopt=wopt) * sin.slp # xRaster - Yrast <- terra::math(sa$aspect, fun="cos", wopt=wopt) * sin.slp # yRaster - Zrast <- terra::math(sa$slope, fun="cos", wopt=wopt) # zRaster - x.sum <- terra::focal(Xrast, w = w, fun=sum, na.rm=na.rm, wopt=wopt) - y.sum <- terra::focal(Yrast, w = w, fun=sum, na.rm=na.rm, wopt=wopt) - z.sum <- terra::focal(Zrast, w = w, fun=sum, na.rm=na.rm, wopt=wopt) - vrm.fun <- function(x, y, z) { - sqrt( (x^2) + (y^2) + (z^2) ) + slope.sin.cos <- terra::lapp(sa$slope, fun = function(s) c(sin(s), cos(s)), wopt = wopt) + aspect.sin.cos <- terra::lapp(sa$aspect, fun = function(a) c(sin(a), cos(a)), wopt = wopt) + + # Decompose into X, Y, Z components + components <- terra::lapp(c(aspect.sin.cos, slope.sin.cos), fun = function(sin_a, cos_a, sin_s, cos_s) { + x <- sin_a * sin_s + y <- cos_a * sin_s + z <- cos_s + c(x, y, z) + }, wopt = wopt) + Xrast <- components[[1]] + Yrast <- components[[2]] + Zrast <- components[[3]] + + # Focal sum for each component + x.sum <- terra::focal(Xrast, w = w, fun = sum, na.rm = na.rm, wopt = wopt) + y.sum <- terra::focal(Yrast, w = w, fun = sum, na.rm = na.rm, wopt = wopt) + z.sum <- terra::focal(Zrast, w = w, fun = sum, na.rm = na.rm, wopt = wopt) + + # Resultant vector + res_vect <- terra::lapp(c(x.sum, y.sum, z.sum), fun = function(x, y, z) sqrt(x^2 + y^2 + z^2), wopt = wopt) + + # Adjust scale factor depending on na.rm + scale.factor <- if (!na.rm) { + round(w[1] * w[2], 0) + } else { + terra::focalCpp(Zrast, w = w, fun = C_CountVals, wopt = wopt) } - res_vect <- terra::lapp(c(x.sum, y.sum, z.sum), fun=vrm.fun, wopt=wopt) #resultant vector - if(!na.rm){ - scale.factor <- round(w[1] * w[2], 0) #Constant scale factor - } else{ - scale.factor<- terra::focalCpp(Zrast, w= w, fun = C_CountVals, wopt=wopt) - } #If include na.rm option, adjust scale.factor to be number of non-NA cells in focal window - out<- 1 - (res_vect / scale.factor) + + # Final VRM output + out <- 1 - (res_vect / scale.factor) names(out)<- "vrm" if(include_scale){names(out)<- paste0(names(out), "_", w[1],"x", w[2])} #Add scale to layer names diff --git a/README.md b/README.md index 9fcdf38..4c5a5fd 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ README ================ Alexander Ilich -March 06, 2025 +April 25, 2025 # MultiscaleDTM @@ -82,6 +82,11 @@ $$ Figure adapted from Walbridge et al., (2018) +- `Pfit` - Calculates multiscale slope and aspect aspect by fitting a + planar surface to the focal window using ordinary least squares. This + will provide equivalent results to using a quadratic fit (Jones, 1998) + but is less computationally expensive. + ### Roughness - `VRM` - Vector ruggedness measure (Sappington et al. 2007) quantifies @@ -407,6 +412,10 @@ Jenness, J.S., 2004. Calculating landscape surface area from digital elevation models. Wildlife Society Bulletin 32, 829–839. +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. + + 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. diff --git a/man/AdjSD.Rd b/man/AdjSD.Rd index ba113d7..3902b34 100644 --- a/man/AdjSD.Rd +++ b/man/AdjSD.Rd @@ -36,9 +36,7 @@ a SpatRaster or RasterLayer of adjusted rugosity 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 (Ilich et al., 2023) } \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) } diff --git a/man/BPI.Rd b/man/BPI.Rd index b4d45c5..54a9569 100644 --- a/man/BPI.Rd +++ b/man/BPI.Rd @@ -42,9 +42,7 @@ A SpatRaster or RasterLayer. Calculates Bathymetric Position Index (BPI). BPI is a measure of relative position that calculates the difference between the value of the focal cell and the mean of cells contained within an annulus shaped neighborhood. Positive values indicate local highs (i.e. peaks) and negative values indicate local lows (i.e. depressions). BPI can be expressed in units of the input DTM raster or can standardized relative to the local topography by dividing by the standard deviation or range of included elevation values in the focal window. BPI calls the function RelPos internally which serves as a general purpose and more flexible function for calculating relative position. } \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) } diff --git a/man/DMV.Rd b/man/DMV.Rd index 8a9d310..2f601f6 100644 --- a/man/DMV.Rd +++ b/man/DMV.Rd @@ -47,9 +47,7 @@ a SpatRaster or RasterLayer. Calculates Difference from Mean Value (DMV). DMV is a measure of relative position that calculates the difference between the value of the focal cell and the mean of all cells in a rectangular or circular neighborhood. Positive values indicate local highs (i.e. peaks) and negative values indicate local lows (i.e. depressions). DMV can be expressed in units of the input DTM raster or can standardized relative to the local topography by dividing by the standard deviation or range of elevation values in the focal window. DMV calls the function RelPos internally which serves as a general purpose and more flexible function for calculating relative position. } \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) } diff --git a/man/DirSlp.Rd b/man/DirSlp.Rd index 19a5411..2899248 100644 --- a/man/DirSlp.Rd +++ b/man/DirSlp.Rd @@ -42,14 +42,16 @@ a SpatRaster or RasterStack of slope and/or aspect (and components of aspect) Calculates the slope along a specified direction. Upslope values are positive and downslope values are negative. } \details{ -dz.dx can be calculated at a specified scale via \code{SlpAsp}, \code{Qfit} (zx and zy), or from an existing layer calculated by another program. +dz.dx and dz.dy can be calculated at a specified scale via \code{SlpAsp}, \code{Pfit}, \code{Qfit} (zx and zy), or from an existing layer calculated by another program. } \examples{ 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) +dirslp3<- DirSlp(alpha = 45, dz.dx= dz3$dz.dx, dz.dy= dz3$dz.dy) } \references{ Neteler, M., & Mitasova, H. (2008). Open source GIS: A GRASS GIS approach (3rd ed.). Springer. diff --git a/man/Pfit.Rd b/man/Pfit.Rd index f9e1c06..1cf658d 100644 --- a/man/Pfit.Rd +++ b/man/Pfit.Rd @@ -2,14 +2,16 @@ % Please edit documentation in R/Pfit.R \name{Pfit} \alias{Pfit} -\title{Calculates multiscale slope and aspect using a local planar fit (IN DEVELOPMENT: DO NOT USE)} +\title{Calculates multiscale slope and aspect using a local planar fit.} \usage{ Pfit( 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() @@ -22,10 +24,14 @@ Pfit( \item{unit}{"degrees" or "radians".} -\item{na.rm}{Logical indicating whether or not to remove NA values before calculations.#'} +\item{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.} + +\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)} +\item{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\emph{pi/2 radians ((-pi/2)- atan2(0,0)+2}pi) and northness and eastness will be calculated from this.} + \item{filename}{character Output filename. Can be a single filename, or as many filenames as there are layers to write a file for each layer} \item{overwrite}{logical. If TRUE, filename is overwritten (default is FALSE).} @@ -36,8 +42,18 @@ Pfit( a SpatRaster (terra) or RasterStack/RasterLayer (raster) } \description{ -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). +} +\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. +} +\examples{ +r<- erupt() +pmetrics<- Pfit(r, w = c(5,5), unit = "degrees", na.rm = TRUE) +plot(pmetrics) } \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. } diff --git a/man/Qfit.Rd b/man/Qfit.Rd index 2e53859..5d20b51 100644 --- a/man/Qfit.Rd +++ b/man/Qfit.Rd @@ -67,9 +67,7 @@ Calculates multiscale slope, aspect, curvature, and morphometric features of a D 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 measured in units of 1/length (e.g. m^-1) except twisting curvature which is measured in radians/length (i.e. change in angle per unit distance), and we adopt a geographic sign convention where convex is positive and concave is negative (i.e., hills are considered convex with positive. 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 (Tg)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. } \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) diff --git a/man/RIE.Rd b/man/RIE.Rd index 5efc7f3..47865dc 100644 --- a/man/RIE.Rd +++ b/man/RIE.Rd @@ -39,9 +39,7 @@ Calculates Roughness Index-Elevation. This is the standard deviation of residual Note the original paper by Cavalli et al (2008) uses a fixed 5x5 window and uses 25 as the denominator indicating use of the population standard deviation. This implementation provides a flexible window size and istead calculates the sample standard deviation which uses a denominator of n-1. } \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) } diff --git a/man/RelPos.Rd b/man/RelPos.Rd index b3b7eea..c9728cb 100644 --- a/man/RelPos.Rd +++ b/man/RelPos.Rd @@ -53,9 +53,7 @@ A SpatRaster or RasterLayer. Calculates the relative position of a focal cell, which represents whether an area is a local high or low. Relative position is the value of the focal cell minus the value of a reference elevation (often the mean of included values in the focal window but see "fun" argument). Positive values indicate local highs (i.e. peaks) and negative values indicate local lows (i.e. depressions). Relative Position can be expressed in units of the input DTM raster or can standardized relative to the local topography by dividing by the standard deviation or range of included elevation values in the focal window. } \examples{ -r<- rast(volcano, extent= ext(2667400, 2667400 + -ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -crs = "EPSG:27200") +r<- erupt() rpos<- RelPos(r, w = c(5,5), shape= "rectangle", exclude_center = TRUE, na.rm = TRUE) plot(rpos) } diff --git a/man/SAPA.Rd b/man/SAPA.Rd index 7dd3103..91074ae 100644 --- a/man/SAPA.Rd +++ b/man/SAPA.Rd @@ -54,9 +54,7 @@ Calculates surface area (Jenness, 2004) to planar area rugosity and by default c Planar area is calculated as the x_dis * y_dis if uncorrected for slope and (x_dis * y_dis)/cos(slope) if corrected for slope. When w=1, this is called "native" scale and is equivalent to what is presented in Du Preez (2015) and available in the ArcGIS Benthic Terrain Modeller add-on. In this case operations are performed on a per cell basis where x_dis is the resolution of the raster in the x direction (left/right) and y_dis is the resolution of the raster in the y direction (up/down) and slope is calculated using the Horn (1981) method. To expand this to multiple scales of analysis, at w > 1 slope is calculated based on Misiuk et al (2021) which provides a modification of the Horn method to extend the matric to multiple spatial scales. Planar area is calculated the same way as for w=1 except that now x_dis is the x resolution of the raster * the number of columns in the focal window, and y_dis is y resolution of the raster * the number of rows. For w > 1, surface area is calculated as the sum of surface areas within the focal window. Although the (modified) Horn slope is used by default to be consistent with Du Preez (2015), slope calculated using a different algorithm (e.g. Wood 1996) could be supplied using the slope_layer argument. Additionally, a slope raster can be supplied if you have already calculated it and do not wish to recalculate it. However, be careful to supply a slope layer measured in radians and calculated at the relevant scale (2 larger than the w of SAPA). } \examples{ -r<- rast(volcano, extent= ext(2667400, 2667400 + -ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -crs = "EPSG:27200") +r<- erupt() sapa<- SAPA(r, w=c(5,5), slope_correction = TRUE) plot(sapa) } diff --git a/man/SlpAsp.Rd b/man/SlpAsp.Rd index 3a1d6c1..ca9a1b7 100644 --- a/man/SlpAsp.Rd +++ b/man/SlpAsp.Rd @@ -54,9 +54,7 @@ dz.dy is the difference in between the weighted mean of the top side of the foca The cells used in these computations is dependent on the "method" chosen. For methods "queen" and "boundary", corner cells have half the weight of all other cells used in the computations. } \examples{ -r<- rast(volcano, extent= ext(2667400, 2667400 + -ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -crs = "EPSG:27200") +r<- erupt() slp_asp<- SlpAsp(r = r, w = c(5,5), unit = "degrees", method = "queen", metrics = c("slope", "aspect", "eastness", "northness")) diff --git a/man/SurfaceArea.Rd b/man/SurfaceArea.Rd index b82d5d7..384028a 100644 --- a/man/SurfaceArea.Rd +++ b/man/SurfaceArea.Rd @@ -30,9 +30,7 @@ a SpatRaster or RasterLayer Calculates surface area on a per cell basis of a DTM based on Jenness, 2004. } \examples{ -r<- rast(volcano, extent= ext(2667400, 2667400 + -ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -crs = "EPSG:27200") +r<- erupt() sa<- SurfaceArea(r) plot(sa) } diff --git a/man/TPI.Rd b/man/TPI.Rd index 2b4f7cf..1a949a7 100644 --- a/man/TPI.Rd +++ b/man/TPI.Rd @@ -47,9 +47,7 @@ SpatRaster or RasterLayer. Calculates Topographic Position Index (TPI). TPI is a measure of relative position that calculates the the difference between the value of the focal cell and the mean of mean of the surrounding cells (i.e. local mean but excluding the value of the focal cell).Positive values indicate local highs (i.e. peaks) and negative values indicate local lows (i.e. depressions). TPI can be expressed in units of the input DTM raster or can standardized relative to the local topography by dividing by the standard deviation or range of included elevation values in the focal window. } \examples{ -r<- rast(volcano, extent= ext(2667400, 2667400 + -ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -crs = "EPSG:27200") +r<- erupt() tpi<- TPI(r, w=c(5,5), shape="rectangle", stand="none", na.rm = TRUE) plot(tpi) } diff --git a/man/VRM.Rd b/man/VRM.Rd index f9b65d7..362bf3d 100644 --- a/man/VRM.Rd +++ b/man/VRM.Rd @@ -39,9 +39,7 @@ Implementation of the Sappington et al., (2007) vector ruggedness measure, modif If the crs is cartesian, when na.rm=TRUE, NA's will be removed from the slope/aspect calculations. When the crs is lat/lon, na.rm=TRUE will not affect the calculation of slope/aspect as terra::terrain will be used since it can calculate slope and aspect for spherical geometry but it does not support na.rm. In both cases when na.rm=TRUE, the x, y, and z components will be summed with na.rm=TRUE, and the N used in the denominator of the VRM equation will be the number of non-NA cells in the window rather than the total number of cells. } \examples{ -r<- rast(volcano, extent= ext(2667400, 2667400 + -ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), -crs = "EPSG:27200") +r<- erupt() vrm<- VRM(r, w=c(5,5), na.rm = TRUE) plot(vrm) } diff --git a/man/fragments/README_Frag.Rmd b/man/fragments/README_Frag.Rmd index 5058032..b29fc64 100644 --- a/man/fragments/README_Frag.Rmd +++ b/man/fragments/README_Frag.Rmd @@ -43,6 +43,8 @@ $$ Figure adapted from Walbridge et al., (2018) +- `Pfit` - Calculates multiscale slope and aspect aspect by fitting a planar surface to the focal window using ordinary least squares. This will provide equivalent results to using a quadratic fit (Jones, 1998) but is less computationally expensive. + ### Roughness - `VRM` - Vector ruggedness measure (Sappington et al. 2007) quantifies roughness by measuring the dispersion of vectors normal to the terrain surface. This is accomplished by calculating the local (3 x 3 cell) slope and aspect, and constructing unit vectors normal to each cell in the DTM. These unit vectors are then decomposed into their corresponding x, y, and z components (i.e. the x, y, and z coordinates of the head of the vector relative to its origin) and used in the following equation (note: N is the number of cells in the window). VRM ranges from zero to one, representing completely smooth to rough surfaces, respectively. @@ -354,6 +356,8 @@ Ilich, A. R., Misiuk, B., Lecours, V., & Murawski, S. A. (2023). MultiscaleDTM: 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 +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 + 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 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. https://doi.org/10.1080/01490410600738021 diff --git a/src/AdjSD.cpp b/src/AdjSD.cpp new file mode 100644 index 0000000..0dccab9 --- /dev/null +++ b/src/AdjSD.cpp @@ -0,0 +1,55 @@ +#define ARMA_WARN_LEVEL 1 +#include +// [[Rcpp::depends(RcppArmadillo)]] +using namespace Rcpp; +using namespace arma; + +//SD from planar fit with na.rm=TRUE +// [[Rcpp::export]] +arma::vec C_AdjSD_narmT(const arma::vec& z, const arma::mat& X_full, int ni, int nw) { + arma::vec out(ni, arma::fill::value(NA_REAL)); + unsigned int thresh = 4; + + for (int i = 0; i < ni; ++i) { + int start = i * nw; + arma::vec zw_full = z.subvec(start, start + nw - 1); + arma::uvec non_na_idx = arma::find_finite(zw_full); + + if (non_na_idx.n_elem >= thresh) { + arma::vec zw = zw_full.elem(non_na_idx); + + arma::vec unique_vals = arma::unique(zw); + if (unique_vals.n_elem == 1) { + out[i] = 0; + } else { + arma::mat X = X_full.rows(non_na_idx); + arma::vec resid = zw - X * arma::solve(X, zw); + out[i] = arma::stddev(resid); + } + } + } + return out; +} + +//SD from planar fit with na.rm=FALSE +// [[Rcpp::export]] +arma::vec C_AdjSD_narmF(const arma::vec& z, const arma::mat& X, const arma::mat& Xt, const arma::mat& XtX_inv, int ni, int nw) { + arma::vec out(ni, arma::fill::value(NA_REAL)); + + for (int i = 0; i < ni; ++i) { + int start = i * nw; + arma::vec Z = z.subvec(start, start + nw - 1); + + if (!Z.has_nan()) { + arma::vec unique_vals = arma::unique(Z); + if (unique_vals.n_elem == 1) { + out[i] = 0; + } else { + arma::vec resid = Z - X * (XtX_inv * (Xt * Z)); + out[i] = arma::stddev(resid); + } + } + } + + return out; +} \ No newline at end of file diff --git a/src/CountVals.cpp b/src/CountVals.cpp new file mode 100644 index 0000000..c2a0077 --- /dev/null +++ b/src/CountVals.cpp @@ -0,0 +1,23 @@ +#include +using namespace Rcpp; + +// [[Rcpp::export]] +NumericVector C_CountVals(NumericVector z, size_t ni, size_t nw) { + NumericVector out(ni, NA_REAL); + const double* z_ptr = z.begin(); // Direct pointer access to z + + for (size_t i = 0; i < ni; ++i) { + size_t start = i * nw; + size_t count = 0; + for (size_t j = 0; j < nw; ++j) { + double val = z_ptr[start + j]; + if (!std::isnan(val)) { + ++count; + } + } + + out[i] = count; + } + + return out; +} diff --git a/src/Qfit.cpp b/src/Qfit.cpp new file mode 100644 index 0000000..c528434 --- /dev/null +++ b/src/Qfit.cpp @@ -0,0 +1,227 @@ +#define ARMA_WARN_LEVEL 1 +#include +// [[Rcpp::depends(RcppArmadillo)]] +using namespace Rcpp; +using namespace arma; + + +//Fit Wood/Evans Quadratic Surface or Plane with Intercept + +//na.rm=TRUE, force_center=FALSE +// [[Rcpp::export]] +arma::mat C_Qfit1_narmT(const arma::vec& z, + const arma::mat& X_full, + bool return_intercept, + int ni, int nw) { + + unsigned int thresh = X_full.n_cols; + + int nlyr = return_intercept ? X_full.n_cols : X_full.n_cols - 1; + + arma::mat out(ni, nlyr, arma::fill::value(NA_REAL)); + + for (int i = 0; i < ni; ++i) { + int start = i * nw; + arma::vec zw = z.subvec(start, start + nw - 1); + arma::uvec valid_idx = arma::find_finite(zw); + + if (valid_idx.n_elem >= thresh) { + arma::vec zw_clean = zw.elem(valid_idx); + arma::mat X = X_full.rows(valid_idx); + + arma::vec unique_vals = arma::unique(zw_clean); + if (unique_vals.n_elem == 1) { + out.row(i).zeros(); + if (return_intercept) { + out(i, nlyr - 1) = unique_vals[0]; // Last coefficient (e.g., intercept "f") + } + } else { + // Solve using least squares: B = solve(X, Z) + arma::vec coef = arma::solve(X, zw_clean); + if (return_intercept) { + out.row(i) = coef.t(); // Full coefficients (including intercept) + } else { + out.row(i) = coef.subvec(0, coef.n_elem - 2).t(); // Exclude intercept + }} + } + } + return out; +} + +//na.rm=FALSE, force_center=FALSE +// [[Rcpp::export]] +arma::mat C_Qfit1_narmF(const arma::vec& z, + const arma::mat& X, + const arma::mat& Xt, + const arma::mat& XtX_inv, + bool return_intercept, + int ni, int nw) { + + int nlyr = return_intercept ? X.n_cols : X.n_cols - 1; + arma::mat out(ni, nlyr, arma::fill::value(NA_REAL)); + + for (int i = 0; i < ni; ++i) { + int start = i * nw; + arma::vec zw = z.subvec(start, start + nw - 1); + + if (!zw.has_nan()) { + arma::vec unique_vals = arma::unique(zw); + + if (unique_vals.n_elem == 1) { + out.row(i).zeros(); + if (return_intercept) { + out(i, nlyr - 1) = unique_vals[0]; // Last coefficient (e.g., intercept "f") + } + } else { + arma::vec coef = XtX_inv * (Xt * zw); + if (return_intercept) { + out.row(i) = coef.t(); // Full coefficients (including intercept) + } else { + out.row(i) = coef.subvec(0, coef.n_elem - 2).t(); // Exclude intercept + }} + } + } + + return out; +} + + +//Fit Wood/Evans Quadratic Surface or Plane forced through center + +//na.rm=TRUE, force_center=TRUE +// [[Rcpp::export]] +arma::mat C_Qfit2_narmT(const arma::vec& z, + const arma::mat& X_full, + int ni, int nw) { + + int nlyr = X_full.n_cols; + arma::mat out(ni, nlyr, arma::fill::value(NA_REAL)); + + unsigned int thresh = nlyr; //Setting to nlayer makes it work for planar fit too + + for (int i = 0; i < ni; ++i) { + int start = i * nw; + arma::vec zw = z.subvec(start, start + nw - 1); + + // Centering the window values + double center_val = zw(nw / 2); + zw -= center_val; + + arma::uvec valid_idx = arma::find_finite(zw); + if (valid_idx.n_elem >= thresh) { + arma::vec zw_clean = zw.elem(valid_idx); + arma::mat X = X_full.rows(valid_idx); + + arma::vec unique_vals = arma::unique(zw_clean); + if (unique_vals.n_elem == 1) { + out.row(i).zeros(); + } else { + // Least squares solution + arma::vec coef = arma::solve(X, zw_clean); + out.row(i) = coef.t(); // row vector + } + } + } + + return out; +} + +//na.rm=FALSE, force_center=TRUE +// [[Rcpp::export]] +arma::mat C_Qfit2_narmF(const arma::vec& z, + const arma::mat& X, + const arma::mat& Xt, + const arma::mat& XtX_inv, + int ni, int nw) { + + int nlyr = X.n_cols; + arma::mat out(ni, nlyr, arma::fill::value(NA_REAL)); + + for (int i = 0; i < ni; ++i) { + int start = i * nw; + arma::vec Z = z.subvec(start, start + nw - 1); + + // Center Z using middle value + double center_val = Z(Z.n_elem / 2); + Z -= center_val; + + if (!Z.has_nan()) { + arma::vec unique_vals = arma::unique(Z); + + if (unique_vals.n_elem == 1) { + out.row(i).zeros(); // All coefficients = 0 + } else { + // OLS coefficients using precomputed XtX_inv and Xt + arma::vec coef = XtX_inv * (Xt * Z); + out.row(i) = coef.t(); + } + } + } + + return out; +} + + +//Planar fit (Not needed. Generalized Qfit) + +//na.rm=TRUE, force_center=FALSE +// arma::mat C_Pfit1_narmT(const arma::vec& z, +// const arma::mat& X_full, +// int ni, int nw) { +// +// int nlyr = X_full.n_cols-1; +// arma::mat out(ni, nlyr, arma::fill::value(NA_REAL)); +// +// unsigned int thresh = 3; +// +// for (int i = 0; i < ni; ++i) { +// int start = i * nw; +// arma::vec zw = z.subvec(start, start + nw - 1); +// arma::uvec valid_idx = arma::find_finite(zw); +// +// if (valid_idx.n_elem >= thresh) { +// arma::vec zw_clean = zw.elem(valid_idx); +// arma::mat X = X_full.rows(valid_idx); +// +// arma::vec unique_vals = arma::unique(zw_clean); +// if (unique_vals.n_elem == 1) { +// out.row(i).zeros(); +// } else { +// // Solve using least squares: B = solve(X, Z) +// arma::vec coef = arma::solve(X, zw_clean); +// out.row(i) = coef.subvec(0, coef.n_elem - 2).t(); // Transpose to write as row +// } +// } +// } +// +// return out; +// } + +//na.rm=FALSE, force_center=FALSE (Not Used) +// arma::mat C_Pfit1_narmF(const arma::vec& z, +// const arma::mat& X, +// const arma::mat& Xt, +// const arma::mat& XtX_inv, +// int ni, int nw) { +// +// int nlyr = X.n_cols - 1; +// arma::mat out(ni, nlyr, arma::fill::value(NA_REAL)); +// +// for (int i = 0; i < ni; ++i) { +// int start = i * nw; +// arma::vec zw = z.subvec(start, start + nw - 1); +// +// if (!zw.has_nan()) { +// arma::vec unique_vals = arma::unique(zw); +// +// if (unique_vals.n_elem == 1) { +// out.row(i).zeros(); +// } else { +// arma::vec coef = XtX_inv * (Xt * zw); +// out.row(i) = coef.subvec(0, coef.n_elem - 2).t(); // Transpose to write as row +// } +// } +// } +// return out; +// } + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 570cdf3..08fd698 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,176 +11,111 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif -// subset_mat_rows -Rcpp::NumericMatrix subset_mat_rows(Rcpp::NumericMatrix Input_Matrix, Rcpp::LogicalVector Input_Log_Vec); -RcppExport SEXP _MultiscaleDTM_subset_mat_rows(SEXP Input_MatrixSEXP, SEXP Input_Log_VecSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Input_Matrix(Input_MatrixSEXP); - Rcpp::traits::input_parameter< Rcpp::LogicalVector >::type Input_Log_Vec(Input_Log_VecSEXP); - rcpp_result_gen = Rcpp::wrap(subset_mat_rows(Input_Matrix, Input_Log_Vec)); - return rcpp_result_gen; -END_RCPP -} -// C_OLS_params -NumericVector C_OLS_params(arma::mat X, arma::mat Y); -RcppExport SEXP _MultiscaleDTM_C_OLS_params(SEXP XSEXP, SEXP YSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Y(YSEXP); - rcpp_result_gen = Rcpp::wrap(C_OLS_params(X, Y)); - return rcpp_result_gen; -END_RCPP -} -// C_OLS_params2 -NumericVector C_OLS_params2(arma::mat Xt, arma::mat XtX_inv, arma::mat Y); -RcppExport SEXP _MultiscaleDTM_C_OLS_params2(SEXP XtSEXP, SEXP XtX_invSEXP, SEXP YSEXP) { +// C_AdjSD_narmT +arma::vec C_AdjSD_narmT(const arma::vec& z, const arma::mat& X_full, int ni, int nw); +RcppExport SEXP _MultiscaleDTM_C_AdjSD_narmT(SEXP zSEXP, SEXP X_fullSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type Xt(XtSEXP); - Rcpp::traits::input_parameter< arma::mat >::type XtX_inv(XtX_invSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Y(YSEXP); - rcpp_result_gen = Rcpp::wrap(C_OLS_params2(Xt, XtX_inv, Y)); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X_full(X_fullSEXP); + Rcpp::traits::input_parameter< int >::type ni(niSEXP); + Rcpp::traits::input_parameter< int >::type nw(nwSEXP); + rcpp_result_gen = Rcpp::wrap(C_AdjSD_narmT(z, X_full, ni, nw)); return rcpp_result_gen; END_RCPP } -// C_OLS_resid -NumericVector C_OLS_resid(arma::mat X, arma::mat Y); -RcppExport SEXP _MultiscaleDTM_C_OLS_resid(SEXP XSEXP, SEXP YSEXP) { +// C_AdjSD_narmF +arma::vec C_AdjSD_narmF(const arma::vec& z, const arma::mat& X, const arma::mat& Xt, const arma::mat& XtX_inv, int ni, int nw); +RcppExport SEXP _MultiscaleDTM_C_AdjSD_narmF(SEXP zSEXP, SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Y(YSEXP); - rcpp_result_gen = Rcpp::wrap(C_OLS_resid(X, Y)); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Xt(XtSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type XtX_inv(XtX_invSEXP); + Rcpp::traits::input_parameter< int >::type ni(niSEXP); + Rcpp::traits::input_parameter< int >::type nw(nwSEXP); + rcpp_result_gen = Rcpp::wrap(C_AdjSD_narmF(z, X, Xt, XtX_inv, ni, nw)); return rcpp_result_gen; END_RCPP } -// C_OLS_resid2 -NumericVector C_OLS_resid2(arma::mat X, arma::mat Xt, arma::mat XtX_inv, arma::mat Y); -RcppExport SEXP _MultiscaleDTM_C_OLS_resid2(SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP YSEXP) { +// C_CountVals +NumericVector C_CountVals(NumericVector z, size_t ni, size_t nw); +RcppExport SEXP _MultiscaleDTM_C_CountVals(SEXP zSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xt(XtSEXP); - Rcpp::traits::input_parameter< arma::mat >::type XtX_inv(XtX_invSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Y(YSEXP); - rcpp_result_gen = Rcpp::wrap(C_OLS_resid2(X, Xt, XtX_inv, Y)); + Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); + Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); + Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); + rcpp_result_gen = Rcpp::wrap(C_CountVals(z, ni, nw)); return rcpp_result_gen; END_RCPP } // C_Qfit1_narmT -NumericMatrix C_Qfit1_narmT(NumericVector z, NumericMatrix X_full, size_t ni, size_t nw); -RcppExport SEXP _MultiscaleDTM_C_Qfit1_narmT(SEXP zSEXP, SEXP X_fullSEXP, SEXP niSEXP, SEXP nwSEXP) { +arma::mat C_Qfit1_narmT(const arma::vec& z, const arma::mat& X_full, bool return_intercept, int ni, int nw); +RcppExport SEXP _MultiscaleDTM_C_Qfit1_narmT(SEXP zSEXP, SEXP X_fullSEXP, SEXP return_interceptSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type X_full(X_fullSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_Qfit1_narmT(z, X_full, ni, nw)); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X_full(X_fullSEXP); + Rcpp::traits::input_parameter< bool >::type return_intercept(return_interceptSEXP); + Rcpp::traits::input_parameter< int >::type ni(niSEXP); + Rcpp::traits::input_parameter< int >::type nw(nwSEXP); + rcpp_result_gen = Rcpp::wrap(C_Qfit1_narmT(z, X_full, return_intercept, ni, nw)); return rcpp_result_gen; END_RCPP } // C_Qfit1_narmF -NumericMatrix C_Qfit1_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, size_t ni, size_t nw); -RcppExport SEXP _MultiscaleDTM_C_Qfit1_narmF(SEXP zSEXP, SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP niSEXP, SEXP nwSEXP) { +arma::mat C_Qfit1_narmF(const arma::vec& z, const arma::mat& X, const arma::mat& Xt, const arma::mat& XtX_inv, bool return_intercept, int ni, int nw); +RcppExport SEXP _MultiscaleDTM_C_Qfit1_narmF(SEXP zSEXP, SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP return_interceptSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xt(XtSEXP); - Rcpp::traits::input_parameter< arma::mat >::type XtX_inv(XtX_invSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_Qfit1_narmF(z, X, Xt, XtX_inv, ni, nw)); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Xt(XtSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type XtX_inv(XtX_invSEXP); + Rcpp::traits::input_parameter< bool >::type return_intercept(return_interceptSEXP); + Rcpp::traits::input_parameter< int >::type ni(niSEXP); + Rcpp::traits::input_parameter< int >::type nw(nwSEXP); + rcpp_result_gen = Rcpp::wrap(C_Qfit1_narmF(z, X, Xt, XtX_inv, return_intercept, ni, nw)); return rcpp_result_gen; END_RCPP } // C_Qfit2_narmT -NumericMatrix C_Qfit2_narmT(NumericVector z, NumericMatrix X_full, size_t ni, size_t nw); +arma::mat C_Qfit2_narmT(const arma::vec& z, const arma::mat& X_full, int ni, int nw); RcppExport SEXP _MultiscaleDTM_C_Qfit2_narmT(SEXP zSEXP, SEXP X_fullSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type X_full(X_fullSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X_full(X_fullSEXP); + Rcpp::traits::input_parameter< int >::type ni(niSEXP); + Rcpp::traits::input_parameter< int >::type nw(nwSEXP); rcpp_result_gen = Rcpp::wrap(C_Qfit2_narmT(z, X_full, ni, nw)); return rcpp_result_gen; END_RCPP } // C_Qfit2_narmF -NumericMatrix C_Qfit2_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, size_t ni, size_t nw); +arma::mat C_Qfit2_narmF(const arma::vec& z, const arma::mat& X, const arma::mat& Xt, const arma::mat& XtX_inv, int ni, int nw); RcppExport SEXP _MultiscaleDTM_C_Qfit2_narmF(SEXP zSEXP, SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xt(XtSEXP); - Rcpp::traits::input_parameter< arma::mat >::type XtX_inv(XtX_invSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type z(zSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Xt(XtSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type XtX_inv(XtX_invSEXP); + Rcpp::traits::input_parameter< int >::type ni(niSEXP); + Rcpp::traits::input_parameter< int >::type nw(nwSEXP); rcpp_result_gen = Rcpp::wrap(C_Qfit2_narmF(z, X, Xt, XtX_inv, ni, nw)); return rcpp_result_gen; END_RCPP } -// C_AdjSD_narmT -NumericVector C_AdjSD_narmT(NumericVector z, NumericMatrix X_full, size_t ni, size_t nw); -RcppExport SEXP _MultiscaleDTM_C_AdjSD_narmT(SEXP zSEXP, SEXP X_fullSEXP, SEXP niSEXP, SEXP nwSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type X_full(X_fullSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_AdjSD_narmT(z, X_full, ni, nw)); - return rcpp_result_gen; -END_RCPP -} -// C_AdjSD_narmF -NumericVector C_AdjSD_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, size_t ni, size_t nw); -RcppExport SEXP _MultiscaleDTM_C_AdjSD_narmF(SEXP zSEXP, SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP niSEXP, SEXP nwSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xt(XtSEXP); - Rcpp::traits::input_parameter< arma::mat >::type XtX_inv(XtX_invSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_AdjSD_narmF(z, X, Xt, XtX_inv, ni, nw)); - return rcpp_result_gen; -END_RCPP -} -// C_Pfit1_narmF -NumericMatrix C_Pfit1_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, LogicalVector idx, size_t ni, size_t nw); -RcppExport SEXP _MultiscaleDTM_C_Pfit1_narmF(SEXP zSEXP, SEXP XSEXP, SEXP XtSEXP, SEXP XtX_invSEXP, SEXP idxSEXP, SEXP niSEXP, SEXP nwSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xt(XtSEXP); - Rcpp::traits::input_parameter< arma::mat >::type XtX_inv(XtX_invSEXP); - Rcpp::traits::input_parameter< LogicalVector >::type idx(idxSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_Pfit1_narmF(z, X, Xt, XtX_inv, idx, ni, nw)); - return rcpp_result_gen; -END_RCPP -} // C_TriArea double C_TriArea(double a, double b, double c); RcppExport SEXP _MultiscaleDTM_C_TriArea(SEXP aSEXP, SEXP bSEXP, SEXP cSEXP) { @@ -195,12 +130,12 @@ BEGIN_RCPP END_RCPP } // C_SurfaceArea -NumericVector C_SurfaceArea(NumericVector z, double x_res, double y_res, bool na_rm, size_t ni, size_t nw); +NumericVector C_SurfaceArea(const NumericVector& z, double x_res, double y_res, bool na_rm, size_t ni, size_t nw); RcppExport SEXP _MultiscaleDTM_C_SurfaceArea(SEXP zSEXP, SEXP x_resSEXP, SEXP y_resSEXP, SEXP na_rmSEXP, SEXP niSEXP, SEXP nwSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); + Rcpp::traits::input_parameter< const NumericVector& >::type z(zSEXP); Rcpp::traits::input_parameter< double >::type x_res(x_resSEXP); Rcpp::traits::input_parameter< double >::type y_res(y_resSEXP); Rcpp::traits::input_parameter< bool >::type na_rm(na_rmSEXP); @@ -210,36 +145,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// C_CountVals -NumericVector C_CountVals(NumericVector z, size_t ni, size_t nw); -RcppExport SEXP _MultiscaleDTM_C_CountVals(SEXP zSEXP, SEXP niSEXP, SEXP nwSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type z(zSEXP); - Rcpp::traits::input_parameter< size_t >::type ni(niSEXP); - Rcpp::traits::input_parameter< size_t >::type nw(nwSEXP); - rcpp_result_gen = Rcpp::wrap(C_CountVals(z, ni, nw)); - return rcpp_result_gen; -END_RCPP -} static const R_CallMethodDef CallEntries[] = { - {"_MultiscaleDTM_subset_mat_rows", (DL_FUNC) &_MultiscaleDTM_subset_mat_rows, 2}, - {"_MultiscaleDTM_C_OLS_params", (DL_FUNC) &_MultiscaleDTM_C_OLS_params, 2}, - {"_MultiscaleDTM_C_OLS_params2", (DL_FUNC) &_MultiscaleDTM_C_OLS_params2, 3}, - {"_MultiscaleDTM_C_OLS_resid", (DL_FUNC) &_MultiscaleDTM_C_OLS_resid, 2}, - {"_MultiscaleDTM_C_OLS_resid2", (DL_FUNC) &_MultiscaleDTM_C_OLS_resid2, 4}, - {"_MultiscaleDTM_C_Qfit1_narmT", (DL_FUNC) &_MultiscaleDTM_C_Qfit1_narmT, 4}, - {"_MultiscaleDTM_C_Qfit1_narmF", (DL_FUNC) &_MultiscaleDTM_C_Qfit1_narmF, 6}, - {"_MultiscaleDTM_C_Qfit2_narmT", (DL_FUNC) &_MultiscaleDTM_C_Qfit2_narmT, 4}, - {"_MultiscaleDTM_C_Qfit2_narmF", (DL_FUNC) &_MultiscaleDTM_C_Qfit2_narmF, 6}, {"_MultiscaleDTM_C_AdjSD_narmT", (DL_FUNC) &_MultiscaleDTM_C_AdjSD_narmT, 4}, {"_MultiscaleDTM_C_AdjSD_narmF", (DL_FUNC) &_MultiscaleDTM_C_AdjSD_narmF, 6}, - {"_MultiscaleDTM_C_Pfit1_narmF", (DL_FUNC) &_MultiscaleDTM_C_Pfit1_narmF, 7}, + {"_MultiscaleDTM_C_CountVals", (DL_FUNC) &_MultiscaleDTM_C_CountVals, 3}, + {"_MultiscaleDTM_C_Qfit1_narmT", (DL_FUNC) &_MultiscaleDTM_C_Qfit1_narmT, 5}, + {"_MultiscaleDTM_C_Qfit1_narmF", (DL_FUNC) &_MultiscaleDTM_C_Qfit1_narmF, 7}, + {"_MultiscaleDTM_C_Qfit2_narmT", (DL_FUNC) &_MultiscaleDTM_C_Qfit2_narmT, 4}, + {"_MultiscaleDTM_C_Qfit2_narmF", (DL_FUNC) &_MultiscaleDTM_C_Qfit2_narmF, 6}, {"_MultiscaleDTM_C_TriArea", (DL_FUNC) &_MultiscaleDTM_C_TriArea, 3}, {"_MultiscaleDTM_C_SurfaceArea", (DL_FUNC) &_MultiscaleDTM_C_SurfaceArea, 6}, - {"_MultiscaleDTM_C_CountVals", (DL_FUNC) &_MultiscaleDTM_C_CountVals, 3}, {NULL, NULL, 0} }; diff --git a/src/SurfaceArea.cpp b/src/SurfaceArea.cpp new file mode 100644 index 0000000..28e1533 --- /dev/null +++ b/src/SurfaceArea.cpp @@ -0,0 +1,91 @@ +#include +using namespace Rcpp; + + +// Calculate area of triangle based on side lengths +// [[Rcpp::export]] +double C_TriArea (double a, double b, double c){ + double s = (a+b+c)/2; + double out =sqrt(s*(s-a)*(s-b)*(s-c)); + return out; +} + +// [[Rcpp::export]] +NumericVector C_SurfaceArea(const NumericVector& z, double x_res, double y_res, bool na_rm, size_t ni, size_t nw) { + NumericVector out(ni, NA_REAL); + + double Lx2 = x_res * x_res; + double Ly2 = y_res * y_res; + double Ld2 = Lx2 + Ly2; + + const double* z_ptr = z.begin(); + + for (size_t i = 0; i < ni; ++i) { + const double* zw = z_ptr + i * nw; + + // Fast access to elevation values + double A = zw[0], B = zw[1], C = zw[2], D = zw[3], E = zw[4], + F = zw[5], G = zw[6], H = zw[7], I = zw[8]; + //|A|B|C| + //|D|E|F| + //|G|H|I| + + // Horizontal + double AB = sqrt(Lx2 + pow(A - B, 2)) / 2; + double BC = sqrt(Lx2 + pow(B - C, 2)) / 2; + double DE = sqrt(Lx2 + pow(D - E, 2)) / 2; + double EF = sqrt(Lx2 + pow(E - F, 2)) / 2; + double GH = sqrt(Lx2 + pow(G - H, 2)) / 2; + double HI = sqrt(Lx2 + pow(H - I, 2)) / 2; + // Vertical + double AD = sqrt(Ly2 + pow(A - D, 2)) / 2; + double BE = sqrt(Ly2 + pow(B - E, 2)) / 2; + double CF = sqrt(Ly2 + pow(C - F, 2)) / 2; + double DG = sqrt(Ly2 + pow(D - G, 2)) / 2; + double EH = sqrt(Ly2 + pow(E - H, 2)) / 2; + double FI = sqrt(Ly2 + pow(F - I, 2)) / 2; + + // Diagonal + double EA = sqrt(Ld2 + pow(E - A, 2)) / 2; + double EC = sqrt(Ld2 + pow(E - C, 2)) / 2; + double EG = sqrt(Ld2 + pow(E - G, 2)) / 2; + double EI = sqrt(Ld2 + pow(E - I, 2)) / 2; + + // Triangles + double tri[8] = { + C_TriArea(EA, AB, BE), + C_TriArea(BE, BC, EC), + C_TriArea(AD, DE, EA), + C_TriArea(EC, CF, EF), + C_TriArea(DE, DG, EG), + C_TriArea(EF, FI, EI), + C_TriArea(EG, EH, GH), + C_TriArea(EH, EI, HI) + }; + + if (na_rm) { + double sum = 0; + int count = 0; + for (int t = 0; t < 8; ++t) { + if (!std::isnan(tri[t])) { + sum += tri[t]; + count++; + } + } + if (count > 0) out[i] = (sum / count) * 8; + } else { + bool has_na = false; + for (int t = 0; t < 8; ++t) { + if (std::isnan(tri[t])) { + has_na = true; + break; + } + } + if (!has_na) { + out[i] = tri[0] + tri[1] + tri[2] + tri[3] + tri[4] + tri[5] + tri[6] + tri[7]; + } + } + } + + return out; +} \ No newline at end of file diff --git a/src/cpp_code.cpp b/src/cpp_code.cpp deleted file mode 100644 index 5458bb1..0000000 --- a/src/cpp_code.cpp +++ /dev/null @@ -1,347 +0,0 @@ -#define ARMA_WARN_LEVEL 1 -#include -// [[Rcpp::depends(RcppArmadillo)]] -using namespace Rcpp; -using namespace arma; - - -//Subsets rows in a numeric matrix according to a logical vector -//https://stackoverflow.com/questions/59284212/efficient-matrix-subsetting-with-rcpp -// [[Rcpp::export]] -Rcpp::NumericMatrix subset_mat_rows(Rcpp::NumericMatrix Input_Matrix, - Rcpp::LogicalVector Input_Log_Vec) { - - // Get the number of rows and columns of "Input_Matrix" - int Mat_nrows = Input_Matrix.nrow(); - int Mat_ncols = Input_Matrix.ncol(); - - // Define matrix "Output_Mat" with the aimed for dimensions, - // i.e. with the number of rows corresponding to the length of "Input_Log_Vec" - Rcpp::NumericMatrix Output_Mat(sum(Input_Log_Vec), Mat_ncols); - - // Loop over the entries of "Input_Log_Vec" to see whether the - // corresponding row of "Input_Matrix" should be included in - // the output matrix "Output_Mat" - for (int i = 0, j = 0; i < Mat_nrows; i++) { - if (Input_Log_Vec[i]) { - Output_Mat(j, _) = Input_Matrix(i, _); - j = j+1; - } - } - - // Return the output matrix "Output_Mat" - return(Output_Mat); -} - -//Ordinary Least Squares (only returns parameters) -// [[Rcpp::export]] -NumericVector C_OLS_params(arma::mat X, arma::mat Y){ - NumericVector B = Rcpp::as(wrap(solve(X,Y))); - return B; - } - -// [[Rcpp::export]] -NumericVector C_OLS_params2(arma::mat Xt, arma::mat XtX_inv, arma::mat Y){ - NumericVector B = Rcpp::as(wrap(XtX_inv * (Xt * Y))); - return B; - } - -//Ordinary Least Squares (only returns residuals) -// [[Rcpp::export]] -NumericVector C_OLS_resid(arma::mat X, arma::mat Y){ - arma::mat B = solve(X,Y); - arma::mat Yhat = X*B; - NumericVector resid = Rcpp::as(wrap(Yhat - Y)); - return resid; - } - -// [[Rcpp::export]] -NumericVector C_OLS_resid2(arma::mat X, arma::mat Xt, arma::mat XtX_inv, arma::mat Y){ - arma::mat B = XtX_inv * (Xt * Y); - arma::mat Yhat = X*B; - NumericVector resid = Rcpp::as(wrap(Yhat - Y)); - return resid; -} - -//Fit Wood/Evans Quadratic Surface with Intercept - -//na.rm=TRUE, force_center=FALSE -// [[Rcpp::export]] -NumericMatrix C_Qfit1_narmT(NumericVector z, NumericMatrix X_full, size_t ni, size_t nw) { - - size_t nlyr = X_full.ncol(); //Number of layers - NumericMatrix out(ni, nlyr); - out.fill(NA_REAL); - colnames(out)= CharacterVector::create("a", "b", "c", "d", "e", "f"); - - int thresh = 6; //NEED AT LEAST 6 POINTS TO CALCULATE BECAUSE NEED AS MANY POINTS AS PARAMETERS - - for (size_t i=0; i(X), as(Z)); - } - }} - return out; -} - -//na.rm=FALSE, force_center=FALSE -// [[Rcpp::export]] -NumericMatrix C_Qfit1_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, size_t ni, size_t nw) { - - size_t nlyr = X.n_cols; //Number of layers - NumericMatrix out(ni, nlyr); - out.fill(NA_REAL); - colnames(out)= CharacterVector::create("a", "b", "c", "d", "e", "f"); - - for (size_t i=0; i::epsilon(); - size_t nlyr = X_full.ncol(); //Number of layers - NumericMatrix out(ni, nlyr); - out.fill(NA_REAL); - colnames(out)= CharacterVector::create("a", "b", "c", "d", "e"); - - int thresh = 5; //NEED AT LEAST 5 POINTS TO CALCULATE BECAUSE NEED AS MANY POINTS AS PARAMETERS - - for (size_t i=0; i(X), as(Z)); - } - }} - return out; -} - -//na.rm=FALSE, force_center=TRUE -// [[Rcpp::export]] -NumericMatrix C_Qfit2_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, size_t ni, size_t nw) { - size_t nlyr = X.n_cols; //Number of layers - NumericMatrix out(ni, nlyr); - out.fill(NA_REAL); - colnames(out)= CharacterVector::create("a", "b", "c", "d", "e"); - - for (size_t i=0; i(X), as(Z)); - out[i] = sd(resid); //SD resid - } - }} - return out; -} - -//SD from planar fit with na.rm=FALSE -// [[Rcpp::export]] -NumericVector C_AdjSD_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, size_t ni, size_t nw){ - NumericVector out(ni, NA_REAL); - //Z = dX + eY + f - - for (size_t i=0; i< ni; i++) { - size_t start = i*nw; - size_t end = start+nw-1; - arma::mat Z = z[Rcpp::Range(start,end)]; //Current window of elevation values - if(Z.has_nan()) {} else { - arma::vec uni_Zvals = unique(Z); - if(uni_Zvals.size() == 1){ - out[i] = 0; //If all Z values are the same residuals are 0 - } else{ - NumericVector resid = C_OLS_resid2(X, Xt, XtX_inv, Z); - out[i] = sd(resid); //SD resid - } - }} - return out; -} - -//Planar fit with na.rm=FALSE -// [[Rcpp::export]] -NumericMatrix C_Pfit1_narmF(NumericVector z, arma::mat X, arma::mat Xt, arma::mat XtX_inv, LogicalVector idx, size_t ni, size_t nw) { - - size_t nlyr = 3; //Number of layers - NumericMatrix out(ni, nlyr); - out.fill(NA_REAL); - colnames(out)= CharacterVector::create("d", "e", "f"); - - for (size_t i=0; i values(mat=TRUE) |> `dimnames<-`(NULL) + expect<- Qfit(r, w = c(5,7), unit = "degrees", metrics = c("qslope", "qaspect", "qeastness", "qnorthness"), na.rm = FALSE, outlier_quantile = c(0,1)) |> values(mat=TRUE) |> `dimnames<-`(NULL) + expect_equal(test, expect) +}) + test_that("Test VRM", { test<- VRM(r, w=c(5,7), na.rm = TRUE) |> values(mat=TRUE) expect<- readRDS(system.file("testdata", "vrm.RDS", package= "MultiscaleDTM"))