Merge pull request #485 from OldLipe/hotfix/bbox

Hotfix in `sits_regularize` function

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: sits
 Type: Package
-Version: 0.15.1
+Version: 0.15.1-1
 Title: Satellite Image Time Series Analysis for Remote Sensing Data Cubes
 Authors@R: c(person('Rolf', 'Simoes', role = c('aut'), email = 'rolf.simoes@inpe.br'),
              person('Gilberto', 'Camara', role = c('aut', 'cre'), email = 'gilberto.camara@inpe.br'),

NAMESPACE

@@ -11,6 +11,8 @@ S3method(.raster_crs,raster)
 S3method(.raster_crs,terra)
 S3method(.raster_data_type,raster)
 S3method(.raster_data_type,terra)
+S3method(.raster_extent,raster)
+S3method(.raster_extent,terra)
 S3method(.raster_extract,raster)
 S3method(.raster_extract,terra)
 S3method(.raster_focal,raster)
@@ -83,6 +85,7 @@ S3method(.source_items_new,bdc_cube)
 S3method(.source_items_new,deafrica_cube)
 S3method(.source_items_new,usgs_cube)
 S3method(.source_items_new,wtss_cube)
+S3method(.source_items_tile_get_bbox,aws_cube)
 S3method(.source_items_tile_get_bbox,stac_cube)
 S3method(.source_items_tile_get_crs,"mspc_cube_landsat-8-c2-l2")
 S3method(.source_items_tile_get_crs,"mspc_cube_sentinel-2-l2a")

NEWS.md

@@ -4,6 +4,10 @@ We are preparing to release the package on CRAN and are making relevant changes
 
 # What's new in SITS version 0.15.1
 
+### New features in SITS version 0.15.1-1
+* Fix bug in `sits_regularize` that generated images with artifacts
+* Fix wrong *bbox* in `sits_cube` from STAC AWS Sentinel-2
+
 ### New features in SITS version 0.15.1
 * Update README.Rmd
 * Support `sits_timeline()` to sits model objects

R/sits_config.R

@@ -977,6 +977,21 @@ sits_list_collections <- function(source = NULL) {
     return(res)
 }
 
+#' @title Retrieve the gdalcubes chunk size
+#' @name .config_gdalcubes_chunk_size
+#' @keywords internal
+#' @return a numeric vector with chunk size
+.config_gdalcubes_chunk_size <- function() {
+
+    chunk_size <- .config_get(key = c("gdalcubes_chunk_size"))
+
+    # post-condition
+    .check_num(chunk_size, len_min = 3, len_max = 3,
+               msg = "invalid 'gdalcubes_chunk_size' in config file")
+
+    return(chunk_size)
+}
+
 #' @title Retrieve the raster package to be used
 #' @name .config_raster_pkg
 #' @keywords internal

R/sits_gdalcubes.R

@@ -12,18 +12,19 @@
 #' @param output_dir  Directory where the aggregated images will be written.
 #' @param cloud_mask  A \code{logical} corresponds to the use of the cloud band
 #'  for aggregation.
+#' @param multicores  A \code{numeric} with the number of cores will be used in
+#'  the regularize. By default is used 1 core.
 #' @param ...         Additional parameters that can be included. See
 #'  '?gdalcubes::write_tif'.
-#' @param version     A \code{character} with version of the output files.
 #'
 #' @return  A data cube tile with information used in its creation.
 .gc_new_cube <- function(tile,
                          cv,
                          img_col,
                          path_db,
                          output_dir,
-                         cloud_mask, ...,
-                         version = "v1") {
+                         cloud_mask,
+                         multicores, ...) {
 
     # set caller to show in errors
     .check_set_caller(".gc_new_cube")
@@ -61,13 +62,43 @@
             collection = .cube_collection(cube = tile)
         )
 
-        list(type   = format_type,
-             nodata = .cube_band_missing_value(cube = cube, band = band),
-             scale  = 1,
-             offset = 0
+        return(
+            list(type   = format_type,
+                 nodata = .cube_band_missing_value(cube = cube, band = band),
+                 scale  = 1,
+                 offset = 0
+            )
         )
     }
 
+    .get_cube_chunks <- function(cv) {
+
+        bbox <- c(xmin = cv[["space"]][["left"]],
+                  xmax = cv[["space"]][["right"]],
+                  ymin = cv[["space"]][["bottom"]],
+                  ymax = cv[["space"]][["top"]])
+
+        size_x <- (max(bbox[c("xmin", "xmax")]) - min(bbox[c("xmin", "xmax")]))
+        size_y <- (max(bbox[c("ymin", "ymax")]) - min(bbox[c("ymin", "ymax")]))
+
+        # a vector with time, x and y
+        chunk_size <- .config_gdalcubes_chunk_size()
+
+        chunks_x <- round(size_x / cv[["space"]][["dx"]]) / chunk_size[[2]]
+        chunks_y <- round(size_y / cv[["space"]][["dy"]]) / chunk_size[[3]]
+
+        # guaranteeing that it will return fewer blocks than calculated
+        num_chunks <- (ceiling(chunks_x) * ceiling(chunks_y)) - 1
+
+        return(max(1, num_chunks))
+    }
+
+    # setting threads to process
+    # multicores number must be smaller than chunks
+    gdalcubes::gdalcubes_options(
+        threads = min(multicores, .get_cube_chunks(cv))
+    )
+
     for (band in .cube_bands(tile, add_cloud = FALSE)) {
 
         # create a raster_cube object to each band the select below change
@@ -149,22 +180,16 @@
 #'  about the cube brick metadata.
 .gc_raster_cube <- function(cube, img_col, cv, cloud_mask) {
 
-
-    chunk_size <- .config_get(key = "gdalcubes_chunk_size")
-
-
     mask_band <- NULL
     if (cloud_mask)
         mask_band <- .gc_cloud_mask(cube)
 
-    gdalcubes_chunk_size <- .config_get(key = "gdalcubes_chunk_size")
-
     # create a brick of raster_cube object
     cube_brick <- gdalcubes::raster_cube(
         image_collection = img_col,
         view = cv,
         mask = mask_band,
-        chunking = gdalcubes_chunk_size)
+        chunking = .config_gdalcubes_chunk_size())
 
     return(cube_brick)
 }
@@ -369,13 +394,13 @@
 
     # create a list of cube view
     cv <- gdalcubes::cube_view(
-        extent = list(left   = roi$left,
-                      right  = roi$right,
-                      bottom = roi$bottom,
-                      top    = roi$top,
+        extent = list(left   = roi[["left"]],
+                      right  = roi[["right"]],
+                      bottom = roi[["bottom"]],
+                      top    = roi[["top"]],
                       t0 = format(toi[["max_min_date"]], "%Y-%m-%d"),
                       t1 = format(toi[["min_max_date"]], "%Y-%m-%d")),
-        srs = tile$crs[[1]],
+        srs = tile[["crs"]][[1]],
         dt = period,
         dx = res,
         dy = res,

R/sits_raster_api.R

@@ -485,6 +485,15 @@
     UseMethod(".raster_crs", pkg_class)
 }
 
+#' @name .raster_properties
+.raster_extent <- function(r_obj, ...) {
+
+    # check package
+    pkg_class <- .raster_check_package()
+
+    UseMethod(".raster_extent", pkg_class)
+}
+
 #' @title Raster package internal frequency values function
 #' @name .raster_freq
 #' @keywords internal

R/sits_raster_api_raster.R

@@ -255,6 +255,15 @@
     )
 }
 
+#' @keywords internal
+#' @export
+.raster_extent.raster <- function(r_obj, ...) {
+
+    suppressWarnings(
+        as.vector(raster::extent(x = r_obj))
+    )
+}
+
 #' @keywords internal
 #' @export
 .raster_freq.raster <- function(r_obj, ...) {

R/sits_raster_api_terra.R

@@ -279,6 +279,15 @@
     )
 }
 
+#' @keywords internal
+#' @export
+.raster_extent.terra <- function(r_obj, ...) {
+
+    suppressWarnings(
+        as.vector(terra::ext(x = r_obj))
+    )
+}
+
 #' @keywords internal
 #' @export
 .raster_freq.terra <- function(r_obj, ...) {

R/sits_regularize.R

@@ -76,8 +76,8 @@
 sits_regularize <- function(cube,
                             output_dir,
                             period,
-                            res        = NULL,
-                            roi        = NULL,
+                            res,
+                            roi = NULL,
                             agg_method = "median",
                             resampling = "bilinear",
                             cloud_mask = TRUE,
@@ -96,51 +96,58 @@ sits_regularize <- function(cube,
                 msg = "sits_regularize not available for collection ",
                 cube$collection, " from ", cube$source
     )
+
     # precondition - test if provided object is a raster cube
     .check_that(
         x = inherits(cube, "raster_cube"),
         msg = paste("provided cube is invalid,",
                     "please provide a 'raster_cube' object.",
                     "see '?sits_cube' for more information.")
     )
+
     # precondition - check output dir fix
     output_dir <- normalizePath(output_dir)
     # verifies the path to save the images
     .check_that(
         x = dir.exists(output_dir),
         msg = "invalid 'output_dir' parameter."
     )
+
     # append gdalcubes path
     path_db <- paste0(output_dir, "/gdalcubes.db")
+
     # precondition - is the period valid?
     .check_na(lubridate::duration(period), msg = "invalid period specified")
 
     # precondition - is the resolution valid?
-    # is there a single resolution? Is not, "res" needs to be set
-    if (length(unique(cube$file_info[[1]]$res)) > 1) {
-        .check_num(x = res,
-                   allow_zero = FALSE,
-                   min = 1,
-                   len_min = 1,
-                   len_max = 1,
-                   msg = "a valid resolution needs to be provided")
-    }
+    .check_num(x = res,
+               allow_zero = FALSE,
+               min = 1,
+               len_min = 1,
+               len_max = 1,
+               msg = "a valid resolution needs to be provided")
+
     # precondition - is the aggregation valid?
     .check_chr_within(
         x = agg_method,
         within = .config_get("gdalcubes_aggreg_methods"),
         discriminator = "any_of",
         msg = "invalid aggregation method"
     )
+
     # precondition - is the resampling valid?
     .check_chr_within(
         x = resampling,
         within = .config_get("gdalcubes_resampling_methods"),
         discriminator = "any_of",
         msg = "invalid resampling method"
     )
+
     # precondition - is the cloud mask valid?
-    .check_lgl_type(cloud_mask, msg = "cloud mask parameter should be TRUE/FALSE")
+    .check_lgl_type(
+        cloud_mask, msg = "cloud mask parameter should be TRUE/FALSE"
+    )
+
     # if the cloud mask is true, is there a cloud band?
     if (cloud_mask) {
         .check_chr_contains(
@@ -149,17 +156,16 @@ sits_regularize <- function(cube,
             msg = "no cloud band available in the cube"
         )
     }
+
     # precondition - is the multicores valid?
     .check_num(
         x = multicores,
-        allow_zero = FALSE,
         min = 1,
+        len_min = 1,
+        len_max = 1,
         msg = "invalid 'multicores' parameter."
     )
 
-    # setting in global env multicores options
-    gdalcubes::gdalcubes_options(threads = multicores)
-
     if (!is.null(roi)) {
 
         # filter only intersecting tiles
@@ -213,13 +219,15 @@ sits_regularize <- function(cube,
                                    agg_method = agg_method,
                                    resampling = resampling)
 
+
         # create of the aggregate cubes
         gc_tile <- .gc_new_cube(tile = tile,
                                 cv = cv,
                                 img_col = img_col,
                                 path_db = path_db,
                                 output_dir = output_dir,
-                                cloud_mask = cloud_mask)
+                                cloud_mask = cloud_mask,
+                                multicores = multicores)
         return(gc_tile)
 
     })

R/sits_source_api_aws.R

@@ -102,3 +102,22 @@
 
     return(crs)
 }
+
+#' @keywords internal
+#' @export
+.source_items_tile_get_bbox.aws_cube <- function(source, ...,
+                                                 tile_items,
+                                                 collection = NULL) {
+    r_obj <- .raster_open_rast(
+        .source_item_get_hrefs(source = source,
+                               item = tile_items$features[[1]])[[1]]
+    )
+
+    # get image bbox
+    bbox <- .raster_extent(r_obj)
+
+    if (is.null(names(bbox)))
+        names(bbox) <- c("xmin", "xmax", "ymin", "ymax")
+
+    return(bbox)
+}