diff --git a/NAMESPACE b/NAMESPACE
index 188225161..ebcfc7907 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -325,7 +325,6 @@ S3method(sits_apply,raster_cube)
 S3method(sits_apply,sits)
 S3method(sits_as_sf,raster_cube)
 S3method(sits_as_sf,sits)
-S3method(sits_bands,base_raster_cube)
 S3method(sits_bands,default)
 S3method(sits_bands,patterns)
 S3method(sits_bands,raster_cube)
diff --git a/R/api_check.R b/R/api_check.R
index 79a6105c8..230109e2f 100644
--- a/R/api_check.R
+++ b/R/api_check.R
@@ -1766,7 +1766,7 @@
     n_bands <- length(.samples_bands.sits(samples))
     n_times <- length(.samples_timeline(samples))
     if(inherits(samples, "sits_base"))
-        n_bands_base <- length(.samples_bands_base(samples))
+        n_bands_base <- length(.samples_base_bands(samples))
     else
         n_bands_base <- 0
     .check_that(ncol(pred) == 2 + n_bands * n_times + n_bands_base)
diff --git a/R/api_classify.R b/R/api_classify.R
index f0c25ab00..520bdadc1 100755
--- a/R/api_classify.R
+++ b/R/api_classify.R
@@ -17,7 +17,9 @@
 #' in the classified images for each corresponding year.
 #'
 #' @param  tile            Single tile of a data cube.
-#' @param  band            Band to be produced.
+#' @param  out_band        Band to be produced.
+#' @param  bands           Bands to extract time series
+#' @param  base_bands      Base bands to extract values
 #' @param  ml_model        Model trained by \code{\link[sits]{sits_train}}.
 #' @param  block           Optimized block to be read into memory.
 #' @param  roi             Region of interest.
@@ -29,7 +31,9 @@
 #' @param  progress        Show progress bar?
 #' @return List of the classified raster layers.
 .classify_tile <- function(tile,
-                           band,
+                           out_band,
+                           bands,
+                           base_bands,
                            ml_model,
                            block,
                            roi,
@@ -42,7 +46,7 @@
     # Output file
     out_file <- .file_derived_name(
         tile = tile,
-        band = band,
+        band = out_band,
         version = version,
         output_dir = output_dir
     )
@@ -53,7 +57,7 @@
         }
         probs_tile <- .tile_derived_from_file(
             file = out_file,
-            band = band,
+            band = out_band,
             base_tile = tile,
             labels = .ml_labels_code(ml_model),
             derived_class = "probs_cube",
@@ -105,7 +109,8 @@
         values <- .classify_data_read(
             tile = tile,
             block = block,
-            bands = .ml_bands(ml_model),
+            bands = bands,
+            base_bands = base_bands,
             ml_model = ml_model,
             impute_fn = impute_fn,
             filter_fn = filter_fn
@@ -138,7 +143,7 @@
         # Prepare probability to be saved
         band_conf <- .conf_derived_band(
             derived_class = "probs_cube",
-            band = band
+            band = out_band
         )
         offset <- .offset(band_conf)
         if (.has(offset) && offset != 0) {
@@ -181,7 +186,7 @@
     # Merge blocks into a new probs_cube tile
     probs_tile <- .tile_derived_merge_blocks(
         file = out_file,
-        band = band,
+        band = out_band,
         labels = .ml_labels_code(ml_model),
         base_tile = tile,
         block_files = block_files,
@@ -374,11 +379,12 @@
 #' @param  tile            Input tile to read data.
 #' @param  block           Bounding box in (col, row, ncols, nrows).
 #' @param  bands           Bands to extract time series
+#' @param  base_bands      Base bands to extract values
 #' @param  ml_model        Model trained by \code{\link[sits]{sits_train}}.
 #' @param  impute_fn       Imputation function
 #' @param  filter_fn       Smoothing filter function to be applied to the data.
 #' @return A matrix with values for classification.
-.classify_data_read <- function(tile, block, bands,
+.classify_data_read <- function(tile, block, bands, base_bands,
                                 ml_model, impute_fn, filter_fn) {
     # For cubes that have a time limit to expire (MPC cubes only)
     tile <- .cube_token_generator(tile)
@@ -388,7 +394,7 @@
         tile = tile,
         block = block
     )
-    # Read and preprocess values of each band
+    # Read and preprocess values of each eo band
     values <- purrr::map(bands, function(band) {
         # Get band values (stops if band not found)
         values <- .tile_read_block(
@@ -436,9 +442,23 @@
         # Return values
         return(as.data.frame(values))
     })
+    # Read and preprocess values of each base band
+    values_base <- purrr::map(base_bands, function(band) {
+        # Read and preprocess values of each base band
+        values_base <- .tile_read_block(
+            tile = .tile_base_info(tile),
+            band = band,
+            block = block
+        )
+        # Return values
+        return(as.data.frame(values_base))
+    })
+    # Combine two lists
+    values <- c(values, values_base)
     # collapse list to get data.frame
-    values <- suppressMessages(purrr::list_cbind(values,
-                                                 name_repair = "universal"))
+    values <- suppressMessages(
+        purrr::list_cbind(values, name_repair = "universal")
+    )
     # Compose final values
     values <- as.matrix(values)
     # Set values features name
diff --git a/R/api_cube.R b/R/api_cube.R
index b08301533..5e61a1483 100644
--- a/R/api_cube.R
+++ b/R/api_cube.R
@@ -131,11 +131,29 @@ NULL
     }
     .cube_set_class(cube)
 }
+#' @title Identity function for data cubes
+#' @keywords internal
+#' @noRd
+#' @name .cube
+#' @param x  cube
+#'
+#' @return data cube object.
 .cube <- function(x) {
     # return the cube
     x
 }
-#' @title Return areas of classes of a class_cue
+#' @title Get base info from a data cube
+#' @keywords internal
+#' @noRd
+#' @name .cube
+#' @param x  cube
+#'
+#' @return data cube from base_info
+.cube_base_info <- function(x) {
+    # return base info data cube
+    dplyr::bind_rows(x[["base_info"]])
+}
+#' @title Return areas of classes of a class_cube
 #' @keywords internal
 #' @noRd
 #' @name .cube_class_areas
@@ -192,7 +210,7 @@ NULL
         class(cube) <- c("raster_cube", class(cube))
         bands <- .cube_bands(cube)
     } else {
-        stop(.conf("messages", "cube_bands"))
+        stop(.conf("messages", ".cube_bands"))
     }
     return(bands)
 }
@@ -203,7 +221,7 @@ NULL
         cube <- tibble::as_tibble(cube)
         bands <- .cube_bands(cube, add_cloud, dissolve)
     } else {
-        stop(.conf("messages", "cube_bands"))
+        stop(.conf("messages", ".cube_bands"))
     }
     return(bands)
 }
@@ -544,6 +562,17 @@ NULL
     }
     return(is_regular)
 }
+
+#' @title Check that cube is a base cube
+#' @name .cube_is_base
+#' @keywords internal
+#' @noRd
+#' @param cube  datacube
+#' @return Called for side effects.
+.cube_is_base <- function(cube) {
+    inherits(cube, "base_raster_cube")
+}
+
 #' @title Find out how many images are in cube during a period
 #' @noRd
 #' @param cube  A data cube.
diff --git a/R/api_data.R b/R/api_data.R
index 3c740ecfc..e11ea5669 100644
--- a/R/api_data.R
+++ b/R/api_data.R
@@ -49,6 +49,10 @@
     } else {
         cld_band <- NULL
     }
+    if (.cube_is_base(cube)) {
+        bands <- setdiff(bands, .cube_bands(.cube_base_info(cube)))
+    }
+
     # define parallelization strategy
     # find block size
     rast <- .raster_open_rast(.tile_path(cube))
@@ -81,9 +85,7 @@
     }
     if (.has(cube[["base_info"]])) {
         # get base info
-        cube_base <- cube[["base_info"]]
-        # bind all base info
-        cube_base <- dplyr::bind_rows(cube_base)
+        cube_base <- .cube_base_info(cube)
         # get bands
         bands_base <- .cube_bands(cube_base)
         # extract data
@@ -97,6 +99,7 @@
         )
         # save base data
         ts_tbl[["base_data"]] <- base_tbl[["time_series"]]
+        # add base class
         class(ts_tbl) <- c("sits_base", class(ts_tbl))
     }
     return(ts_tbl)
diff --git a/R/api_plot_time_series.R b/R/api_plot_time_series.R
index 552f35493..01dc1ae5e 100644
--- a/R/api_plot_time_series.R
+++ b/R/api_plot_time_series.R
@@ -74,7 +74,7 @@
             # how many time series are to be plotted?
             number <- nrow(data2)
             # what are the band names?
-            bands <- .samples_bands(data2)
+            bands <- .samples_bands(data2, include_base = FALSE)
             # what are the reference dates?
             ref_dates <- .samples_timeline(data2)
             # align all time series to the same dates
diff --git a/R/api_plot_vector.R b/R/api_plot_vector.R
index 31d3883e7..9d0924cb3 100644
--- a/R/api_plot_vector.R
+++ b/R/api_plot_vector.R
@@ -87,7 +87,7 @@
     # verifies if stars package is installed
     .check_require_packages("stars")
     # verifies if tmap package is installed
-    .check_require_packages("plot")
+    .check_require_packages("tmap")
     # precondition - check color palette
     .check_palette(palette)
     # revert the palette
diff --git a/R/api_predictors.R b/R/api_predictors.R
index 2a5797f6a..779028b2f 100644
--- a/R/api_predictors.R
+++ b/R/api_predictors.R
@@ -19,11 +19,12 @@
     # Get samples time series
     pred <- .ts(samples)
     # By default get bands as the same of first sample
-    bands <- .samples_bands(samples)
+    bands <- .samples_bands(samples, include_base = FALSE)
     # Preprocess time series
     if (.has(ml_model)) {
         # If a model is informed, get predictors from model bands
-        bands <- .ml_bands(ml_model)
+        bands <- intersect(.ml_bands(ml_model), bands)
+
         # Normalize values for old version model classifiers that
         #   do not normalize values itself
         # Models trained after version 1.2 do this automatically before
@@ -49,7 +50,7 @@
             })
         }
     }
-    # Create predictors...
+    # Create predictors
     pred <- pred[c(.pred_cols, bands)]
     # Add sequence 'index' column grouped by 'sample_id'
     pred <- pred |>
@@ -68,36 +69,22 @@
 }
 #' @export
 .predictors.sits_base <- function(samples, ml_model = NULL) {
-    # Get predictors for time series
     # Prune samples time series
     samples <- .samples_prune(samples)
     # Get samples time series
-    pred <- .ts(samples)
-    # By default get bands as the same of first sample
-    bands <- .samples_bands.sits(samples)
-    # Create predictors...
-    pred <- pred[c(.pred_cols, bands)]
-    # Add sequence 'index' column grouped by 'sample_id'
-    pred <- pred |>
-        dplyr::select("sample_id", "label", dplyr::all_of(bands)) |>
-        dplyr::group_by(.data[["sample_id"]]) |>
-        dplyr::mutate(index = seq_len(dplyr::n())) |>
-        dplyr::ungroup()
-    # Rearrange data to create predictors
-    pred <- tidyr::pivot_wider(
-        data = pred, names_from = "index", values_from = dplyr::all_of(bands),
-        names_prefix = if (length(bands) == 1) bands else "",
-        names_sep = ""
-    )
-    # get predictors for base data
-    base <- dplyr::bind_rows(samples$base_data)
-    base <- base[,-1]
-    # join time series predictors with base data predictors
-    pred <- dplyr::bind_cols(pred, base)
+    pred <- .predictors.sits(samples, ml_model)
+    # Get predictors for base data
+    pred_base <- samples |>
+                 dplyr::rename(
+                     "_" = "time_series", "time_series" = "base_data"
+                 ) |>
+                 .predictors.sits() |>
+                 dplyr::select(-.data[["label"]])
+    # Merge predictors
+    pred <- dplyr::inner_join(pred, pred_base, by = "sample_id")
     # Return predictors
     pred
 }
-
 #' @title Get predictors names with timeline
 #' @keywords internal
 #' @noRd
@@ -113,7 +100,6 @@
         USE.NAMES = FALSE
     ))
 }
-
 #' @title Get features from predictors
 #' @keywords internal
 #' @noRd
diff --git a/R/api_samples.R b/R/api_samples.R
index 87d35f576..232a157a2 100644
--- a/R/api_samples.R
+++ b/R/api_samples.R
@@ -89,23 +89,28 @@
 #' @noRd
 #' @param samples Data.frame with samples
 #' @return Bands for the first sample
-.samples_bands <- function(samples) {
+.samples_bands <- function(samples, ...) {
     # Bands of the first sample governs whole samples data
     UseMethod(".samples_bands", samples)
 }
 #' @export
-.samples_bands.sits <- function(samples) {
+.samples_bands.sits <- function(samples, ...) {
     # Bands of the first sample governs whole samples data
     bands <- setdiff(names(.samples_ts(samples)), "Index")
     return(bands)
 }
 #' @export
-.samples_bands.sits_base <- function(samples) {
+.samples_bands.sits_base <- function(samples, include_base = TRUE) {
     # Bands of the first sample governs whole samples data
-    ts_bands <- .samples_bands.sits(samples)
-    base_bands <- .samples_bands_base(samples)
-    bands <- c(ts_bands, base_bands)
-    return(bands)
+    bands <- .samples_bands.sits(samples)
+
+    if (include_base) {
+        bands <- c(
+            bands, .samples_base_bands(samples)
+        )
+    }
+
+    bands
 }
 #' @export
 .samples_bands.default <- function(samples) {
@@ -113,15 +118,28 @@
     ts_bands <- .samples_bands.sits(samples)
     return(ts_bands)
 }
+#' @title Check if samples is base (has base property)
+#' @noRd
+#' @param samples Data.frame with samples
+#' @return TRUE/FALSE
+.samples_is_base <- function(samples) {
+    inherits(samples, "sits_base")
+}
+#' @title Get samples base data (if available)
+#' @noRd
+#' @param samples Data.frame with samples
+#' @return data.frame with base data.
+.samples_base_data <- function(samples) {
+    samples[["base_data"]]
+}
 #' @title Get bands of base data for samples
 #' @noRd
 #' @param samples Data.frame with samples
 #' @return Bands for the first sample
-.samples_bands_base <- function(samples) {
+.samples_base_bands <- function(samples) {
     # Bands of the first sample governs whole samples data
-    setdiff(names(samples$base_data[[1]]), "Index")
+    setdiff(names(samples[["base_data"]][[1]]), "Index")
 }
-
 #' @title Get timeline of time series samples
 #' @noRd
 #' @param samples Data.frame with samples
@@ -129,7 +147,6 @@
 .samples_timeline <- function(samples) {
     as.Date(samples[["time_series"]][[1]][["Index"]])
 }
-
 #' @title Select bands of time series samples
 #' @noRd
 #' @param samples Data.frame with samples
diff --git a/R/api_tile.R b/R/api_tile.R
index 62d5b1afc..bcd5b7102 100644
--- a/R/api_tile.R
+++ b/R/api_tile.R
@@ -534,8 +534,8 @@ NULL
 #' @export
 .tile_bands.base_raster_cube <- function(tile, add_cloud = TRUE) {
     bands <- .tile_bands.raster_cube(tile, add_cloud)
-    base_bands <- .tile_bands_base(tile)
-    all_bands <- c(bands, base_bands)
+    base_bands <- .tile_bands.raster_cube(.tile_base_info(tile))
+    unique(c(bands, base_bands))
 }
 #' @export
 .tile_bands.default <- function(tile, add_cloud = TRUE) {
@@ -544,14 +544,6 @@ NULL
     bands <- .tile_bands(tile, add_cloud)
     return(bands)
 }
-#' @title Get bands of base data for tile
-#' @noRd
-#' @param samples Data.frame with samples
-#' @return Bands for the first sample
-.tile_bands_base <- function(tile) {
-    # Bands of the first sample governs whole samples data
-    names(tile$base_info[[1]])
-}
 #' @title Set bands in tile file_info.
 #' @rdname .tile_bands
 #' @keywords internal
@@ -1617,6 +1609,7 @@ NULL
     })
     return(cog_sizes)
 }
+
 #' @title  Return base info
 #' @name .tile_base_info
 #' @keywords internal
@@ -1625,9 +1618,6 @@ NULL
 #'
 #' @param  tile       Tile to be plotted
 #' @return            Base info tibble
-#'
-#'
 .tile_base_info <- function(tile) {
-    if (.cube_has_base_info(tile))
-        return(tile[["base_info"]][[1]])
+    return(tile[["base_info"]][[1]])
 }
diff --git a/R/sits_add_base_cube.R b/R/sits_add_base_cube.R
index f08c81903..7e294fd9c 100644
--- a/R/sits_add_base_cube.R
+++ b/R/sits_add_base_cube.R
@@ -55,13 +55,13 @@
 #' }
 #' @export
 #'
-sits_add_base_cube <- function(cube1, cube2){
+sits_add_base_cube <- function(cube1, cube2) {
     .check_set_caller("sits_add_base_cube")
     .check_is_raster_cube(cube1)
     .check_that(.cube_is_regular(cube1))
     .check_that(inherits(cube2, "dem_cube"))
     # pre-condition for merge is having the same tiles
-    .check_that(all(cube1[["tile"]] %in% cube2[["tile"]]))
+    .check_cubes_same_tiles(cube1, cube2)
     # extract tiles
     tiles <- .cube_tiles(cube1)
     # add base info by tile
@@ -70,8 +70,11 @@ sits_add_base_cube <- function(cube1, cube2){
         tile_cube2 <- .cube_filter_tiles(cube2, tile_name)
         # get files from 2nd cube
         fi_cube2 <- .fi(tile_cube2)
+        # get timelines
+        tile_cube1_tl <- .tile_timeline(tile_cube1)
+        tile_cube2_tl <- .tile_timeline(tile_cube2)
         # align timelines
-        fi_cube2[["date"]] <- .fi_min_date(.fi(tile_cube1))
+        fi_cube2[["date"]] <- tile_cube1_tl[1:length(tile_cube2_tl)]
         # update 2nd cube files
         .fi(tile_cube2) <- fi_cube2
         # append cube to base info
diff --git a/R/sits_bands.R b/R/sits_bands.R
index c960c8bd5..66dec3bb5 100644
--- a/R/sits_bands.R
+++ b/R/sits_bands.R
@@ -50,19 +50,6 @@ sits_bands.sits <- function(x) {
 }
 #' @rdname sits_bands
 #' @export
-sits_bands.base_raster_cube <- function(x) {
-    # get time series bands
-    ts_bands <- .tile_bands(x)
-    # get base bands
-    tile_base <- .tile_base_info(x)
-    base_bands <- .tile_bands(tile_base)
-    return(list(
-        "time_series_bands" = ts_bands,
-        "base_bands" = base_bands
-    ))
-}
-#' @rdname sits_bands
-#' @export
 sits_bands.raster_cube <- function(x) {
     # set caller to show in errors
     .check_set_caller("sits_bands")
diff --git a/R/sits_classify.R b/R/sits_classify.R
index 6b7cdd3f8..114eab59d 100644
--- a/R/sits_classify.R
+++ b/R/sits_classify.R
@@ -214,7 +214,6 @@ sits_classify.raster_cube <- function(data,
                                       progress = TRUE) {
     # set caller for error messages
     .check_set_caller("sits_classify_raster")
-    # reduce GPU memory for MPS
     # preconditions
     .check_is_raster_cube(data)
     .check_that(.cube_is_regular(data))
@@ -254,19 +253,34 @@ sits_classify.raster_cube <- function(data,
         .check_filter_fn(filter_fn)
     # Retrieve the samples from the model
     samples <- .ml_samples(ml_model)
-    # Retrieve the bands from the model
-    bands <- .ml_bands(ml_model)
+    # By default, base bands is null.
+    base_bands <- NULL
+    if (.cube_is_base(data)) {
+        # Get base bands
+        base_bands <- intersect(
+            .ml_bands(ml_model), .cube_bands(.cube_base_info(data))
+        )
+    }
+    # get non-base bands
+    bands <- setdiff(.ml_bands(ml_model), base_bands)
     # Do the samples and tile match their timeline length?
     .check_samples_tile_match_timeline(samples = samples, tile = data)
     # Do the samples and tile match their bands?
     .check_samples_tile_match_bands(samples = samples, tile = data)
-    # Check memory and multicores
     # Get block size
     block <- .raster_file_blocksize(.raster_open_rast(.tile_path(data)))
     # Check minimum memory needed to process one block
     job_memsize <- .jobs_memsize(
         job_size = .block_size(block = block, overlap = 0),
-        npaths = length(.tile_paths(data, bands)) + length(.ml_labels(ml_model)),
+        npaths = (
+            length(.tile_paths(data, bands)) +
+            length(.ml_labels(ml_model)) +
+            ifelse(
+                test = .cube_is_base(data),
+                yes = length(.tile_paths(.cube_base_info(data), base_bands)),
+                no = 0
+            )
+        ),
         nbytes = 8,
         proc_bloat = proc_bloat
     )
@@ -305,7 +319,9 @@ sits_classify.raster_cube <- function(data,
         # Classify the data
         probs_tile <- .classify_tile(
             tile = tile,
-            band = "probs",
+            out_band = "probs",
+            bands = bands,
+            base_bands = base_bands,
             ml_model = ml_model,
             block = block,
             roi = roi,
diff --git a/R/sits_get_data.R b/R/sits_get_data.R
index ee0053aa6..6f4e6fd9c 100644
--- a/R/sits_get_data.R
+++ b/R/sits_get_data.R
@@ -270,8 +270,7 @@ sits_get_data.sits <- function(cube,
                                impute_fn = impute_linear(),
                                multicores = 2,
                                progress = FALSE) {
-    if (!.has(bands))
-        bands <- .cube_bands(cube)
+    bands <- .default(bands, .cube_bands(cube))
     # Extract time series from a cube given a data.frame
     data <- .data_get_ts(
         cube       = cube,
diff --git a/man/sits_bands.Rd b/man/sits_bands.Rd
index c64ca66fb..f308039c8 100644
--- a/man/sits_bands.Rd
+++ b/man/sits_bands.Rd
@@ -3,7 +3,6 @@
 \name{sits_bands}
 \alias{sits_bands}
 \alias{sits_bands.sits}
-\alias{sits_bands.base_raster_cube}
 \alias{sits_bands.raster_cube}
 \alias{sits_bands.patterns}
 \alias{sits_bands.sits_model}
@@ -18,8 +17,6 @@ sits_bands(x)
 
 \method{sits_bands}{sits}(x)
 
-\method{sits_bands}{base_raster_cube}(x)
-
 \method{sits_bands}{raster_cube}(x)
 
 \method{sits_bands}{patterns}(x)