Merge pull request #114 from bleutner/dev

Merge publication fixes

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: RStoolbox
 Type: Package
 Title: Remote Sensing Data Analysis
-Version: 1.0.0
-Date: 2024-04-24
+Version: 1.0.1
+Date: 2024-08-26
 Authors@R: c(
     person("Benjamin", "Leutner", role= "aut", email="rstoolboxpackage@gmail.com", comment = c(ORCID = "0000-0002-6893-2002")),
     person("Ned", "Horning", role ="aut", email="horning@amnh.org"),
@@ -34,6 +34,7 @@ Imports:
     magrittr
 Suggests:
     randomForest,
+    lattice,
     kernlab,
     e1071,
     gridExtra,
@@ -42,5 +43,5 @@ Suggests:
 LinkingTo: Rcpp,
     RcppArmadillo
 License: GPL (>=3)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 LazyData: true

R/RcppExports.R

@@ -49,7 +49,7 @@ specSimC <- function(x, em) {
     .Call('_RStoolbox_specSimC', PACKAGE = 'RStoolbox', x, em)
 }
 
-spectralIndicesCpp <- function(x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf) {
-    .Call('_RStoolbox_spectralIndicesCpp', PACKAGE = 'RStoolbox', x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf)
+spectralIndicesCpp <- function(x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf, formulas) {
+    .Call('_RStoolbox_spectralIndicesCpp', PACKAGE = 'RStoolbox', x, indices, redBand, blueBand, greenBand, nirBand, redEdge1Band, redEdge2Band, redEdge3Band, swir1Band, swir2Band, swir3Band, maskLayer, maskValue, L, s, G, C1, C2, Levi, swir2ccc, swir2cdiff, sf, formulas)
 }
 

R/internalFunctions.R

@@ -255,7 +255,84 @@
 #' On package startup
 #' @noRd 
 .onLoad <- function(libname, pkgname){
+
+    L <- 0.5
+    G <- 2.5
+    L_evi <- 1
+    C1 <- 6
+    C2 <- 7.5
+    s <- 1
+    swir2ccc <- NULL
+    swir2coc <- NULL
+    NDVI <- NULL
+
+    .IDXDB <- list(
+        CLG = list(c("Gitelson2003", "Green-band Chlorophyll Index"),
+                function(redEdge3, green) {redEdge3/green - 1}),
+        CLRE = list(c("Gitelson2003", "Red-edge-band Chlorophyll Index"),
+                function(redEdge3, redEdge1) {redEdge3/redEdge1 - 1}),
+        CTVI = list(c("Perry1984", "Corrected Transformed Vegetation Index"),
+                function(red, nir) {(NDVI+.5)/sqrt(abs(NDVI+.5))} ),
+        DVI = list(c("Richardson1977", "Difference Vegetation Index"),
+                function(red, nir) {s*nir-red}),
+        EVI = list(c("Huete1999", "Enhanced Vegetation Index"),
+                function(red, nir, blue) {G * ((nir - red) / (nir + C1 * red - C2 * blue + L_evi))}),
+        EVI2 = list(c("Jiang 2008", "Two-band Enhanced Vegetation Index"),
+                function(red, nir) {G * (nir-red)/(nir + 2.4*red +1)}),
+        GEMI = list(c("Pinty1992", "Global Environmental Monitoring Index"),
+                function(red, nir) {(((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)) * (1 - ((((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5) ) / (nir + red + 0.5)) * 0.25)) - ((red - 0.125) / (1 - red))}),
+        GNDVI = list(c("Gitelson1998", "Green Normalised Difference Vegetation Index" ),
+                function(green, nir) {(nir-green)/(nir+green)}),
+		KNDVI = list(c("Camps-Valls2021", "Kernel Normalised Difference Vegetation Index"),
+				function(red, nir) {tanh(((nir-red)/(nir+red)))^2}),
+        MCARI = list(c("Daughtery2000", "Modified Chlorophyll Absorption Ratio Index" ),
+                function(green, red, redEdge1) {((redEdge1-red)-(redEdge1-green))*(redEdge1/red)}),
+        MNDWI = list(c("Xu2006", "Modified Normalised Difference Water Index"),
+                function(green, swir2) {(green-swir2) / (green+swir2)}),
+        MSAVI = list(c("Qi1994", "Modified Soil Adjusted Vegetation Index" ),
+                function(red, nir) {nir + 0.5 - (0.5 * sqrt((2 * nir + 1)^2 - 8 * (nir - (2 * red))))}),
+        MSAVI2 = list(c("Qi1994", "Modified Soil Adjusted Vegetation Index 2" ),
+                function(red, nir) {(2 * (nir + 1) - sqrt((2 * nir + 1)^2 - 8 * (nir - red))) / 2}),
+        MTCI = list(c("DashAndCurran2004", "MERIS Terrestrial Chlorophyll Index"),
+                function(red, redEdge1, redEdge2) {(redEdge2-redEdge1)/(redEdge1-red)}),
+        NBRI = list(c("Garcia1991", "Normalised Burn Ratio Index"),
+                function(nir, swir3) { (nir - swir3) / (nir + swir3)}),
+        NDREI1 = list(c("GitelsonAndMerzlyak1994", "Normalised Difference Red Edge Index 1"),
+                function(redEdge2, redEdge1) {(redEdge2-redEdge1)/(redEdge2+redEdge1)}),
+        NDREI2 = list(c("Barnes2000", "Normalised Difference Red Edge Index 2"),
+                function(redEdge3, redEdge1) {(redEdge3-redEdge1)/(redEdge3+redEdge1)}),
+        NDVI = list(c("Rouse1974", "Normalised Difference Vegetation Index"),
+                function(red, nir) {(nir-red)/(nir+red)}),
+        NDVIC = list(c("Nemani1993", "Corrected Normalised Difference Vegetation Index" ),
+                function(red, nir, swir2) {(nir-red)/(nir+red)*(1-((swir2 - swir2ccc)/(swir2coc-swir2ccc)))}),
+        NDWI = list(c("McFeeters1996", "Normalised Difference Water Index"),
+                function(green, nir) {(green - nir)/(green + nir)}),
+        NDWI2 = list(c("Gao1996", "Normalised Difference Water Index"),
+                function(nir, swir2) {(nir - swir2)/(nir + swir2)}),
+        NRVI = list(c("Baret1991", "Normalised Ratio Vegetation Index" ),
+                function(red, nir) {(red/nir - 1)/(red/nir + 1)}),
+        REIP = list(c("GuyotAndBarnet1988", "Red Edge Inflection Point"),
+                function(red, redEdge1, redEdge2, redEdge3) {0.705+0.35*((red+redEdge3)/(2-redEdge1))/(redEdge2-redEdge1)}),
+        RVI = list(c("", "Ratio Vegetation Index"),
+                function(red, nir) {red/nir}),
+        SATVI = list(c("Marsett2006", "Soil Adjusted Total Vegetation Index"),
+                function(red, swir2, swir3) {(swir2 - red) / (swir2 + red + L) * (1 + L) - (swir3 / 2)}),
+        SAVI = list(c("Huete1988", "Soil Adjusted Vegetation Index"),
+                function(red, nir) {(nir - red) * (1+L) / (nir + red + L)}),
+        SLAVI = list(c("Lymburger2000", "Specific Leaf Area Vegetation Index"),
+                function(red, nir, swir2) {nir / (red + swir2)}),
+        SR = list(c("Birth1968", "Simple Ratio Vegetation Index"),
+                function(red, nir) {nir / red}),
+        TTVI = list(c("Thiam1997", "Thiam's Transformed Vegetation Index"),
+                function(red, nir) {sqrt(abs((nir-red)/(nir+red) + 0.5))}),
+        TVI = list(c("Deering1975", "Transformed Vegetation Index"),
+                function(red, nir) {sqrt((nir-red)/(nir+red)+0.5)}),
+        WDVI = list(c("Richardson1977", "Weighted Difference Vegetation Index"),
+                function(red, nir) {nir - s * red})
+    )
+
     if(is.null(getOption("RStoolbox.verbose")))  options(RStoolbox.verbose = FALSE)
+    if(is.null(getOption("RStoolbox.idxdb")))  options(RStoolbox.idxdb = .IDXDB)
 }
 
 #' Init verbosity within functions 
@@ -270,6 +347,18 @@
     options(RStoolbox.verbose = verbose)
 }
 
+
+#' Init spectralIndices within the index DB
+#'
+#' will restore global options after function has been called
+#' @param spectralIndices List
+#' @keywords internal
+#' @noRd
+.initIDXdb <- function(idxdb){
+    idxbg <- force(getOption("RStoolbox.idxdb"))
+    do.call("on.exit", list(substitute(options(RStoolbox.idxdb = idxbg))), envir=parent.frame())
+    options(RStoolbox.idxdb = idxdb)
+}
 
 
 #' Clean up on package unload

R/mesma.R

@@ -10,7 +10,7 @@
 #' }
 #' @param iterate Integer. Set maximum iteration per pixel. Processing time could increase the more iterations are made possible. Default is 400.
 #' @param tolerance Numeric. Tolerance limit representing a nearly zero minimal number. Default is 1e-8. 
-#' @param n_models Logical. Only applies if \code{em} contains column \code{class}. Defines how many endmember combinations should be picked. Maximum is the minimum number of endmembers of a class. Defaults to 5.
+#' @param n_models Logical. Only applies if \code{em} contains column \code{class}. Defines how many endmember combinations should be picked. Maximum is the minimum number of endmembers of a class. Defaults to the amount of rows of em.
 #' @param sum_to_one Logical. Defines whether a sum-to-one constraint should be applied so that probabilities of endmember classes sum to one (a constraint not covered by NNLS) to be interpretable as fractions. Defaults to \code{TRUE}. To get actual NNLS results, change to \code{FALSE}.
 #' @param verbose Logical. Prints progress messages during execution.
 #' @param ... further arguments passed to \link[terra]{writeRaster}.
@@ -72,7 +72,7 @@
 #' @importFrom terra which.min rast selectRange nlyr
 #' @export
 #' 
-mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.00000001, n_models = 5, sum_to_one = TRUE, ..., verbose){
+mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.00000001, n_models = NULL, sum_to_one = TRUE, ..., verbose){
   img <- .toTerra(img)
   ## messages
   if(!missing("verbose")) .initVerbose(verbose)
@@ -89,6 +89,12 @@ mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.0000000
   if(anyNA(em)){
     stop("'em' is not allowed to contain NA values. Spectra must be consistent.")
   }
+
+  # Check for n_models type
+  if(!is.null(n_models) && !is.numeric(n_models)){
+    stop("'n_models' needs to be an 'numeric' object.")
+  }
+
 
   method <- toupper(method) 
   meth_avail <- c("NNLS") #available methods
@@ -115,6 +121,11 @@ mesma <- function(img, em, method = "NNLS", iterate = 400, tolerance = 0.0000000
     classes <- unique(em$class)
     # check n_models
     n_em_cl <- min(sapply(classes, function(x) nrow(em[em$class == x,]), USE.NAMES = F))
+
+    if(is.null(n_models)){
+      n_models <- n_em_cl
+    }
+
     if(n_em_cl < n_models){
       warning(paste0("'n_models' cannot be larger than the minimum number of provided endmembers per class. Setting 'n_models' to ", n_em_cl, "."))
       n_models <- n_em_cl

R/rsOpts.R

@@ -3,14 +3,20 @@
 #' shortcut to options(RStoolbox.*)
 #' 
 #' @param verbose Logical. If \code{TRUE} many functions will print status messages about the current processing step. By default verbose mode is disabled.
+#' @param idxdb List. The list conatins the formal calculation of spectral indices. Modify this list to pipe your own spectral index through the internal C++ calculation of RStoolbox.
 #' @export
 #' @return
-#' No return, just a setter for the verbosiness of the RStoolbox package
+#' No return, just a setter for the verbosiness and the index-database of the RStoolbox package. For latter, see the example of Rstoolbox::spectralIndices()
 #' @examples 
 #' rsOpts(verbose=TRUE)
 #'
-rsOpts <- function(verbose){
-    options(RStoolbox.verbose=verbose)
+rsOpts <- function(verbose=NULL, idxdb=NULL){
+    if(!is.null(verbose)){
+        options(RStoolbox.verbose=verbose)
+    }
+    if(!is.null(idxdb)){
+        options(RStoolbox.idxdb=idxdb)
+    }
 }