Merge pull request #1074 from OldLipe/feat/appy_prop

Fix legend `sits_reclassify()` bug

Author: gilbertocamara

R/api_apply.R

@@ -149,7 +149,8 @@
         base_tile = feature,
         block_files = block_files,
         multicores = 1,
-        update_bbox = FALSE
+        update_bbox = FALSE,
+        normalized = normalized
     )
     # Return a feature tile
     band_tile

R/api_conf.R

@@ -937,7 +937,7 @@ NULL
 .conf_eo_band <- function(source, collection, band) {
     # Format band name
     band <- .band_eo(band)
-    # Return a default value if band does not exists in config
+    # does the band exists in cube config?
     if (!.conf_eo_band_exists(source, collection, band)) {
         return(NULL)
     }

R/api_raster.R

@@ -652,6 +652,24 @@
     UseMethod(".raster_freq", pkg_class)
 }
 
+#' @title Raster package internal raster data type
+#' @name .raster_datatype
+#' @keywords internal
+#' @noRd
+#' @author Felipe Carvalho, \email{felipe.carvalho@@inpe.br}
+#'
+#' @param r_obj    raster package object
+#' @param by_layer  A logical value indicating the type of return
+#' @param ...      additional parameters to be passed to raster package
+#'
+#' @return A character value with data type
+.raster_datatype <- function(r_obj, ..., by_layer = TRUE) {
+    # check package
+    pkg_class <- .raster_check_package()
+
+    UseMethod(".raster_datatype", pkg_class)
+}
+
 #' @title Raster package internal summary values function
 #' @name .raster_summary
 #' @keywords internal

R/api_raster_terra.R

@@ -459,6 +459,21 @@
     terra::freq(x = r_obj, bylayer = TRUE)
 }
 
+#' @title Raster package internal raster data type
+#' @name .raster_datatype
+#' @keywords internal
+#' @noRd
+#' @author Felipe Carvalho, \email{felipe.carvalho@@inpe.br}
+#'
+#' @param r_obj     raster package object
+#' @param by_layer  A logical value indicating the type of return
+#' @param ...      additional parameters to be passed to raster package
+#'
+#' @return A character value with data type
+.raster_datatype.terra <- function(r_obj, ..., by_layer = TRUE) {
+    terra::datatype(x = r_obj, bylyr = by_layer)
+}
+
 #' @title Summary values of terra object
 #' @keywords internal
 #' @noRd

R/api_reclassify.R

@@ -194,3 +194,19 @@
     # Return closure
     reclassify_fn
 }
+
+.reclassify_new_labels <- function(cube, rules) {
+    # Get cube labels
+    cube_labels <- .cube_labels(cube, dissolve = FALSE)[[1]]
+    # Get rules new labels
+    new_labels <- setdiff(names(rules), cube_labels)
+    # Does rules has new labels in the composition?
+    if (.has(new_labels) > 0) {
+        # Get the next index
+        next_idx <- max(as.numeric(names(cube_labels))) + 1
+        idx_values <- seq.int(
+            from = next_idx, to = next_idx + length(new_labels) - 1                                    )
+        names(new_labels) <- as.character(idx_values)
+    }
+    return(c(cube_labels, new_labels))
+}

R/api_space_time_operations.R

@@ -129,7 +129,7 @@
 .contains <- function(x, y) {
     as_crs <- sf::st_crs(x)
     y <- sf::st_transform(y, crs = as_crs)
-    apply(sf::st_contains(x, y, sparse = FALSE), 1, any)
+    apply(suppressMessages(sf::st_contains(x, y, sparse = FALSE)), 1, any)
 }
 #' @title Find the closest points.
 #'

R/api_tile.R

@@ -526,11 +526,13 @@ NULL
 #' @name .tile_band_conf
 #' @keywords internal
 #' @noRd
-#' @param tile A tile.
-#' @param band Band character vector.
+#' @param tile   A tile.
+#' @param band   Band character vector.
+#' @param normalized A logical indicating if band is normalized.
+#' @param ...    Additional parameters.
 #'
 #' @return band_conf or band_cloud_conf
-.tile_band_conf <- function(tile, band) {
+.tile_band_conf <- function(tile, band, ...) {
     UseMethod(".tile_band_conf", tile)
 }
 #' @export
@@ -552,13 +554,13 @@ NULL
     return(NULL)
 }
 #' @export
-.tile_band_conf.derived_cube <- function(tile, band) {
+.tile_band_conf.derived_cube <- function(tile, band, ...) {
     .conf_derived_band(
         derived_class = .tile_derived_class(tile), band = band[[1]]
     )
 }
 #' @export
-.tile_band_conf.default <- function(tile, band) {
+.tile_band_conf.default <- function(tile, band, ...) {
     tile <- tibble::as_tibble(tile)
     tile <- .cube_find_class(tile)
     band_conf <- .tile_band_conf(tile, band)

R/sits_apply.R

@@ -52,6 +52,12 @@
 #' inside the kernel window. Central pixel is \code{NA} just only
 #' all pixels in the window are \code{NA}.
 #'
+#' By default, the indexes generated by the \code{sits_apply()} function are
+#' normalized between -1 and 1, scaled by a factor of 0.0001.
+#' Normalized indexes are saved as INT2S (Integer with sign).
+#' If the \code{normalized} parameter is FALSE, no scaling factor will be
+#' applied and the index will be saved as FLT4S (Float with sign).
+#'
 #' @section Summarizing kernel functions:
 #' \itemize{
 #' \item{\code{w_median()}: returns the median of the neighborhood's values.}
@@ -125,6 +131,8 @@ sits_apply.raster_cube <- function(data, ...,
     .check_is_regular(data)
     # Check window size
     .check_window_size(window_size)
+    # Check normalized index
+    .check_lgl(normalized)
     # Check memsize
     .check_memsize(memsize, min = 1, max = 16384)
     # Check multicores

R/sits_reclassify.R

@@ -182,8 +182,8 @@ sits_reclassify.class_cube <- function(cube,
             },
             .msg_error = "mask roi does not intersect cube"
         )
-        # Get output labels
-        labels <- unique(c(.cube_labels(cube), names(rules)))
+        # Get new labels from cube and pre-defined rules from user
+        labels <- .reclassify_new_labels(cube, rules)
         # Classify the data
         class_tile <- .reclassify_tile(
             tile = tile,

man/sits-package.Rd

@@ -1,4 +1,4 @@
-test_that("EVI generation", {
+test_that("Testing normalized index generation", {
     s2_cube <- tryCatch(
         {
             sits_cube(
@@ -26,8 +26,8 @@ test_that("EVI generation", {
         suppressWarnings(dir.create(dir_images))
     }
     unlink(list.files(dir_images,
-        pattern = "\\.tif$",
-        full.names = TRUE
+                      pattern = "\\.tif$",
+                      full.names = TRUE
     ))
 
 
@@ -43,12 +43,12 @@ test_that("EVI generation", {
     )
 
     gc_cube_new <- sits_apply(gc_cube,
-        EVI2 = 2.5 * (B8A - B05) / (B8A + 2.4 * B05 + 1),
-        multicores = 1,
-        output_dir = dir_images
+                              EVI = 2.5 * (B8A - B05) / (B8A + 2.4 * B05 + 1),
+                              multicores = 1,
+                              output_dir = dir_images
     )
 
-    expect_true(all(sits_bands(gc_cube_new) %in% c("EVI2", "B05", "B8A")))
+    expect_true(all(sits_bands(gc_cube_new) %in% c("EVI", "B05", "B8A")))
 
     timeline <- sits_timeline(gc_cube_new)
     start_date <- timeline[1]
@@ -63,7 +63,7 @@ test_that("EVI generation", {
     file_info_b8a <- .fi(gc_cube_new) |> .fi_filter_bands(bands = "B8A")
     b8a_band_1 <- .raster_open_rast(file_info_b8a$path[[1]])
 
-    file_info_evi2 <- .fi(gc_cube_new) |> .fi_filter_bands(bands = "EVI2")
+    file_info_evi2 <- .fi(gc_cube_new) |> .fi_filter_bands(bands = "EVI")
     evi2_band_1 <- .raster_open_rast(file_info_evi2$path[[1]])
 
     b05_100 <- as.numeric(b05_band_1[100] / 10000)
@@ -104,14 +104,87 @@ test_that("EVI generation", {
                                 progress = FALSE)
     evi_tibble_2 <- sits_apply(
         evi_tibble,
-        EVI2_NEW = 2.5 * (B8A - B05) / (B8A + 2.4 * B05 + 1)
+        EVI_NEW = 2.5 * (B8A - B05) / (B8A + 2.4 * B05 + 1)
     )
 
-    values_evi2 <- .tibble_time_series(evi_tibble_2)$EVI2
-    values_evi2_new <- .tibble_time_series(evi_tibble_2)$EVI2_NEW
+    values_evi2 <- .tibble_time_series(evi_tibble_2)$EVI
+    values_evi2_new <- .tibble_time_series(evi_tibble_2)$EVI_NEW
     expect_equal(values_evi2, values_evi2_new, tolerance = 0.001)
 })
 
+test_that("Testing non-normalized index generation", {
+    data_dir <- system.file("extdata/raster/mod13q1", package = "sits")
+    cube <- sits_cube(
+        source = "BDC",
+        collection = "MOD13Q1-6",
+        data_dir = data_dir,
+        progress = FALSE
+    )
+
+
+    dir_images <- paste0(tempdir(), "/images/")
+    if (!dir.exists(dir_images)) {
+        suppressWarnings(dir.create(dir_images))
+    }
+    gc_cube_new <- sits_apply(cube,
+                              XYZ = 1 / NDVI * 0.25,
+                              normalized = FALSE,
+                              multicores = 2,
+                              output_dir = dir_images
+    )
+
+    expect_true(all(sits_bands(gc_cube_new) %in% c("NDVI", "XYZ")))
+
+    file_info_ndvi <- .fi(gc_cube_new) |> .fi_filter_bands(bands = "NDVI")
+    ndvi_band_1 <- .raster_open_rast(file_info_ndvi$path[[1]])
+
+    file_info_xyz <- .fi(gc_cube_new) |> .fi_filter_bands(bands = "XYZ")
+    xyz_band_1 <- .raster_open_rast(file_info_xyz$path[[1]])
+
+    scale_factor <- 10000
+    ndvi_100 <- as.numeric(ndvi_band_1[100] / 10000)
+    xyz_100 <- as.numeric(xyz_band_1[100] / 10000) * scale_factor
+
+    xyz_calc_100 <- 1 / ndvi_100 * 0.25
+    expect_equal(xyz_100, xyz_calc_100, tolerance = 0.001)
+
+    ndvi_150 <- as.numeric(ndvi_band_1[150] / 10000)
+    xyz_150 <- as.numeric(xyz_band_1[150] / 10000) * scale_factor
+
+    xyz_calc_150 <- 1 / ndvi_150 * 0.25
+    expect_equal(xyz_150, xyz_calc_150, tolerance = 0.001)
+
+    bbox_cube <- sits_bbox(gc_cube_new, as_crs = "EPSG:4326")
+    lats <- runif(10, min = bbox_cube[["ymin"]], max = bbox_cube[["ymax"]])
+    longs <- runif(10, min = bbox_cube[["xmin"]], max = bbox_cube[["xmax"]])
+
+    timeline <- sits_timeline(gc_cube_new)
+    start_date <- timeline[1]
+    end_date <- timeline[length(timeline)]
+
+    csv_tb <- purrr::map2_dfr(lats, longs, function(lat, long) {
+        tibble::tibble(
+            longitude = long,
+            latitude = lat,
+            start_date = start_date,
+            end_date = end_date,
+            label = "NoClass"
+        )
+    })
+    csv_file <- paste0(tempdir(), "/csv_gc_cube2.csv")
+    write.csv(csv_tb, file = csv_file)
+
+    xyz_tibble <- sits_get_data(gc_cube_new, csv_file, progress = FALSE)
+    xyz_tibble_2 <- sits_apply(
+        xyz_tibble,
+        XYZ_NEW = 1 / NDVI * 0.25
+    )
+
+    values_xyz2 <- .tibble_time_series(xyz_tibble)$XYZ
+    values_xyz_new <- .tibble_time_series(xyz_tibble_2)$XYZ_NEW
+    expect_equal(values_xyz2, values_xyz_new, tolerance = 0.001)
+})
+
 test_that("Kernel functions", {
     data_dir <- system.file("extdata/raster/mod13q1", package = "sits")
     cube <- sits_cube(
@@ -225,8 +298,8 @@ test_that("Kernel functions", {
     expect_true(max_1 == max_2)
 
     tif_files <- grep("tif",
-        list.files(tempdir(), full.names = TRUE),
-        value = TRUE
+                      list.files(tempdir(), full.names = TRUE),
+                      value = TRUE
     )
 
     success <- file.remove(tif_files)