Merge pull request #1418 from M3nin0/fix/area-calculation

review area calculation functions

Author: gilbertocamara

NAMESPACE

@@ -21,6 +21,9 @@ S3method(.check_samples,default)
 S3method(.check_samples,sits)
 S3method(.cube_adjust_crs,default)
 S3method(.cube_adjust_crs,grd_cube)
+S3method(.cube_area_freq,class_cube)
+S3method(.cube_area_freq,class_vector_cube)
+S3method(.cube_area_freq,default)
 S3method(.cube_as_sf,default)
 S3method(.cube_as_sf,raster_cube)
 S3method(.cube_bands,default)

R/api_cube.R

@@ -367,19 +367,13 @@ NULL
 #'
 #' @return A \code{vector} with the areas of the cube labels.
 .cube_class_areas <- function(cube) {
-    # Get area for each class for each row of the cube
-    freq_lst <- slider::slide(cube, function(tile) {
-        # Get the frequency count and value for each labelled image
-        .tile_area_freq(tile)
-    })
-    # Get a tibble by binding the row (duplicated labels with different counts)
-    freq <- do.call(rbind, freq_lst)
+    # Get cube area / pixel frequency
+    freq <- .cube_area_freq(cube)
     # summarize the counts for each label
     freq <- freq |>
         dplyr::filter(!is.na(class)) |>
         dplyr::group_by(class) |>
         dplyr::summarise(area = sum(.data[["area"]]))
-
     # Area is taken as the sum of pixels
     class_areas <- freq[["area"]]
     # Names of area are the classes
@@ -388,7 +382,40 @@ NULL
     class_areas[is.na(class_areas)] <- 0.0
     class_areas
 }
-
+#' @title Return areas/frequency of classes of a class_cube
+#' @keywords internal
+#' @noRd
+#' @name .cube_area_freq
+#' @param cube       class cube
+#'
+#' @return A \code{vector} with the area and pixel frequency of the cube labels.
+.cube_area_freq <- function(cube) {
+    UseMethod(".cube_area_freq", cube)
+}
+#' @export
+.cube_area_freq.class_cube <- function(cube) {
+    # Get area for each class for each row of the cube
+    slider::slide_dfr(cube, function(tile) {
+        # Get the frequency count and value for each labelled image
+        .tile_area_freq(tile)
+    }) |>
+        dplyr::filter(!is.na(.data[["class"]]))
+}
+#' @export
+.cube_area_freq.class_vector_cube <- function(cube) {
+    # Get area for each class for each row of the cube
+    slider::slide_dfr(cube, function(tile) {
+        # Get the frequency count and value for each labelled image
+        .tile_area_freq(tile)
+    }) |>
+        dplyr::filter(!is.na(.data[["class"]]))
+}
+#' @export
+.cube_area_freq.default <- function(cube) {
+    cube <- tibble::as_tibble(cube)
+    cube <- .cube_find_class(cube)
+    .cube_area_freq(cube)
+}
 #' @title Return bands of a data cube
 #' @keywords internal
 #' @noRd

R/api_message.R

@@ -76,6 +76,15 @@
         )
     }
 }
+#' @title Warning when converting a slow area calculation will be performed
+#' @name .message_warnings_slow_area_calculation
+#' @noRd
+#' @returns Called for side effects
+.message_warnings_slow_area_calculation <- function() {
+    if (.message_warnings()) {
+        warning(.conf("messages", "summary_class_cube_slow_area"), call. = FALSE)
+    }
+}
 #' @title Test if progress bar should be shown
 #' @name .message_progress
 #' @noRd

R/api_raster.R

@@ -883,6 +883,20 @@
     terra::freq(x = rast, bylayer = TRUE)
 }
 
+#' @title Raster package internal area value
+#' @name .raster_area
+#' @keywords internal
+#' @noRd
+#' @author Felipe Carlos, \email{efelipecarlos@@gmail.com}
+#'
+#' @param rast    raster package object to count values
+#' @param ...      additional parameters to be passed to raster package
+#'
+#' @return matrix with layer, value, and area columns
+.raster_area <- function(rast, ...) {
+    terra::expanse(x = rast, ...)
+}
+
 #' @title Raster package internal raster data type
 #' @name .raster_datatype
 #' @keywords internal

R/api_tile.R

@@ -1393,24 +1393,53 @@ NULL
 }
 #' @export
 .tile_area_freq.class_cube <- function(tile) {
-    # Open first raster
-    rast <- .raster_open_rast(.tile_path(tile))
-    # Retrieve the frequency
-    freq <- tibble::as_tibble(.raster_freq(rast))
-    # get labels
-    labels <- .tile_labels(tile)
-    # pixel area
-    # convert the area to hectares
-    # assumption: spatial resolution unit is meters
-    area <- freq[["count"]] * .tile_xres(tile) * .tile_yres(tile) / 10000.0
-    # Include class names
-    freq <- dplyr::mutate(
-        freq,
-        area = area,
-        class = labels[as.character(freq[["value"]])]
+    # get tile crs
+    tile_crs <- .tile_crs(tile)
+    # get tile crs unit (metre or degree)
+    tile_crs_unit <- sf::st_crs(tile_crs)$units_gdal
+    # extract the file path
+    tile_file <- .tile_paths(tile)
+    # read the files with terra
+    rast <- .raster_open_rast(tile_file)
+    # get area by pixels
+    if (tile_crs_unit == "metre") {
+        # get a frequency of values
+        class_areas <- .raster_freq(rast) |>
+            dplyr::select(-.data[["layer"]])
+        # transform to km^2
+        cell_size <- .tile_xres(tile) * .tile_yres(tile)
+        class_areas[["area"]] <- (class_areas[["count"]] * cell_size) / 1000000L
+    } else {
+        # get pixels by class
+        class_count <- .raster_freq(rast)
+        # get area by class in km^2
+        class_areas <- .raster_area(rast = rast, unit = "km", byValue = TRUE)
+        # Merge area and pixel count
+        class_areas <- dplyr::full_join(class_count, class_areas, by = "value") |>
+            dplyr::select(-.data[["layer.x"]], -.data[["layer.y"]])
+    }
+    # change value to character
+    class_areas <- dplyr::mutate(
+        class_areas,
+        value = as.character(.data[["value"]])
     )
-    # Return frequencies
-    freq
+    # create a data.frame with the labels
+    tile_labels <- .tile_labels(tile)
+    df1 <- tibble::tibble(
+        value = names(tile_labels),
+        class = unname(tile_labels)
+    )
+    # join the labels with the areas
+    sum_areas <- dplyr::full_join(df1, class_areas, by = "value")
+    sum_areas <- dplyr::mutate(sum_areas,
+                               area = signif(.data[["area"]], 2L),
+                               .keep = "unused"
+    )
+    # replace na
+    sum_clean <- sum_areas |>
+        tidyr::replace_na(list(count = 0L, area = 0.0))
+
+    sum_clean
 }
 #' @export
 .tile_area_freq.class_vector_cube <- function(tile) {

R/sits_summary.R

@@ -447,49 +447,21 @@ summary.variance_cube <- function(object, ...,
 #' @export
 summary.class_cube <- function(object, ...) {
     .check_set_caller("summary_class_cube")
-    # Extract classes values for each tiles using a sample size
-    classes_areas <- slider::slide(object, function(tile) {
-        # extract the file path
-        tile_file <- .tile_paths(tile)
-        # read the files with terra
-        r <- .raster_open_rast(tile_file)
-        # get a frequency of values
-        class_areas <- .raster_freq(r)
-        # transform to km^2
-        cell_size <- .tile_xres(tile) * .tile_yres(tile)
-        class_areas[["area"]] <- (class_areas[["count"]] * cell_size) / 1000000L
-        # change value to character
-        class_areas <- dplyr::mutate(
-            class_areas,
-            value = as.character(.data[["value"]])
-        )
-        # create a data.frame with the labels
-        tile_labels <- .tile_labels(tile)
-        df1 <- tibble::tibble(
-            value = names(tile_labels),
-            class = unname(tile_labels)
-        )
-        # join the labels with the areas
-        sum_areas <- dplyr::full_join(df1, class_areas, by = "value")
-        sum_areas <- dplyr::mutate(sum_areas,
-            area_km2 = signif(.data[["area"]], 2L),
-            .keep = "unused"
-        )
-        # remove layer information
-        sum_clean <- sum_areas[, -3L] |>
-            tidyr::replace_na(list(layer = 1L, count = 0L, area_km2 = 0.0))
-
-        sum_clean
+    # check if cube has only metre crs
+    cube_has_only_metre_crs <- slider::slide_lgl(object, function(tile) {
+        # get tile crs
+        tile_crs <- .tile_crs(tile)
+        # check unit
+        sf::st_crs(tile_crs)$units_gdal == "metre"
     })
-    # Combine tiles areas
-    classes_areas <- dplyr::bind_rows(classes_areas) |>
-        dplyr::group_by(.data[["value"]], .data[["class"]]) |>
-        dplyr::summarise(
-            count = sum(.data[["count"]]),
-            area_km2 = sum(.data[["area_km2"]]),
-            .groups = "keep"
-        ) |>
-        dplyr::ungroup()
-    # Return classes areas
-    classes_areas
+    # all must be true
+    cube_has_only_metre_crs <- all(cube_has_only_metre_crs)
+    # if there is any non-metre crs, inform user about slow area calculation
+    if (!cube_has_only_metre_crs) {
+        # show warning message or slow area calculation
+        .message_warnings_slow_area_calculation()
+    }
+    # Extract cube class areas
+    .cube_area_freq(object) |>
+        dplyr::rename(area_km2 = "area")
 }

inst/extdata/config_messages.yml

@@ -562,5 +562,6 @@ sits_xgboost: "wrong input parameters - see example in documentation"
 summary_raster_cube: "check that input is regular data cube"
 summary_derived_cube: "check that input is probability data cube"
 summary_class_cube: "check that input is classified data cube"
+summary_class_cube_slow_area: "we detected that your cube uses a non-metre Coordinate Reference System (CRS) - this may use more processing time since additional unit transformations will be required."
 summary_class_cube_area: "some classes have no area: "
 expect_error: "called from expect_error"