R/sits_indexes.R

In e-sensing/sits.data: Satellite Image Time Series Data Handling

#' @title Builds tasseled cap bands
#' @name sits_tasseled_cap
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#' @description  Adds new tasseled cap bands.
#' @param data          Valid sits tibble.
#' @param sensor        sensor.
#' @return A sits tibble with same samples and the new bands.
#' @examples
#' # Retrieve data for time series with label samples in Mato Grosso in Brazil
#' data (samples_mt_6bands)
#' # Generate a new image with the tasseled cap
#' tc.tb <- sits_tasseled_cap(samples_mt_6bands, sensor = "MODIS")
#' @export
sits_tasseled_cap <- function(data, sensor = "MODIS"){
    # backward compatibility
    if ("coverage" %in% names(data))
        data <- .sits_tibble_rename(data)
    bands <- sits_bands(data)

    # verify if bands are available
    # get the bands required to do TCAP for the sensor
    b_coef <- .sits_config_tcap(sensor, "brightness")
    bands_b_coef <- names(b_coef)

    # get the original bands
    orig_bands <- .sits_config_original_bands(sensor)
    # keep only the original bands
    bands <- bands[bands %in% orig_bands]

    # Do we have all the bands?
    if (!all(bands_b_coef %in% bands))
        message("Tasseled cap for ", sensor, " needs bands: ",
                    paste0(bands_b_coef, collapse = ", "),"\n",
                    "Data has less bands: ",
                    paste0(bands, collapse = ", "), "\n",
                    "Result will be an approximation")

    # calculate the tasseled cap indexes

    calc_formula <- function(tc){
        # retrieve the transformation coefs
        coef <- .sits_config_tcap(sensor, tc)
        # calculate the formula
        formula.lst <- purrr::pmap(list(bands), function(b) {
            formula <- paste0(coef[b],"*",paste(b))
            return(formula)
        })
        formula <- stringi::stri_join_list(formula.lst , sep = " ",
                                           collapse = " + ")
        return(formula)
    }

    # calculate tasseled cap bands
    formula <- calc_formula("brightness")
    data    <- sits_mutate_bands(data,  tcb = eval(str2lang(formula)))
    formula <- calc_formula("greeness")
    data    <- sits_mutate_bands(data,  tcg = eval(str2lang(formula)))
    formula <- calc_formula("wetness")
    data    <- sits_mutate_bands(data,  tcw = eval(str2lang(formula)))

    return(data)
}

#' @title Builds soil-adjusted vegetation index
#' @name sits_savi
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#' @description  Adds new tasseled cap bands.
#' @param data       Valid sits tibble.
#' @return A sits tibble with the SAVI index.
#' @examples
#' \donttest{
#' # Retrieve data for time series with label samples in Mato Grosso in Brazil
#' data (samples_mt_6bands)
#' # Generate a new image with the tasseled cap
#' savi.tb <- sits_savi(samples_mt_6bands)
#' }
#' @export
sits_savi <- function(data){
    # backward compatibility
    if ("coverage" %in% names(data))
        data <- .sits_tibble_rename(data)
    bands <- sits_bands(data)
    bands_savi <- c("nir", "red")
    assertthat::assert_that(all(bands_savi %in% (bands)),
        msg = "sits_savi: not enough bands to compute")

    data <- sits_mutate_bands(data, savi = (1.5)*(nir - red)/(nir + red + 0.5))

    return(data)
}

#' @title Builds normalized difference water index
#' @name sits_ndwi
#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
#' @description  Adds new tasseled cap bands.
#' @param data       A valid sits tibble.
#' @return A sits tibble with the SAVI index.
#' @examples
#' \donttest{
#' # Retrieve data for time series with label samples in Mato Grosso in Brazil
#' data(samples_mt_6bands)
#' # Generate a new image with the tasseled cap
#' ndwi.tb <- sits_ndwi(samples_mt_6bands)
#' }
#' @export
sits_ndwi <- function(data){
    # backward compatibility
    if ("coverage" %in% names(data))
        data <- .sits_tibble_rename(data)
    bands <- sits_bands(data)
    bands_ndwi <- c("nir", "mir")
    assertthat::assert_that(all(bands_ndwi %in% (bands)),
                         msg = "sits_ndwi: not enough bands to compute")

    data <- sits_mutate_bands(data, ndwi = (1.5) * (nir - mir)/(nir + mir))

    return(data)
}