diff --git a/pysits/docs/content/impute_linear.md b/pysits/docs/content/impute_linear.md
new file mode 100644
index 0000000..f98f927
--- /dev/null
+++ b/pysits/docs/content/impute_linear.md
@@ -0,0 +1,9 @@
+Replace NA values by linear interpolation
+
+Remove NA by linear interpolation
+
+Args:
+    data (list | pandas.DataFrame): A time series vector or matrix.
+
+Returns:
+    R: A set of filtered time series using the imputation function.
diff --git a/pysits/docs/content/impute_mean.md b/pysits/docs/content/impute_mean.md
new file mode 100644
index 0000000..c798bfe
--- /dev/null
+++ b/pysits/docs/content/impute_mean.md
@@ -0,0 +1,9 @@
+Remove NA using mean
+
+Remove NA using mean
+
+Args:
+    data (list[float] | pandas.DataFrame): A time series or matrix.
+
+Returns:
+    R: A set of filtered time series using the imputation function.
diff --git a/pysits/docs/content/impute_mean_window.md b/pysits/docs/content/impute_mean_window.md
new file mode 100644
index 0000000..f995737
--- /dev/null
+++ b/pysits/docs/content/impute_mean_window.md
@@ -0,0 +1,23 @@
+Remove NA using weighted moving average
+
+Remove NA using weighted moving average
+
+Args:
+    data (list): A time series vector or matrix.
+    k (int): Width of the moving average window. Expands to both sides
+        of the center element e.g. k = 2 means 4 observations (2 left,
+        2 right) are taken into account. If all observations in the
+        current window are NA, the window size is automatically
+        increased until there are at least 2 non-NA values present.
+    weighting (str): The weighting strategy to be used. More details
+        below (default is "simple").
+
+Returns:
+    R: A set of filtered time series using the imputation function.
+
+Notes:
+    The `weighting` parameter defines the weighting strategy used in the
+    moving window. The strategies available are:
+    - `simple` - Simple Moving Average (SMA) (default option)
+    - `linear` - Linear Weighted Moving Average (LWMA)
+    - `exponential` - Exponential Weighted Moving Average (EWMA)
diff --git a/pysits/docs/content/impute_median.md b/pysits/docs/content/impute_median.md
new file mode 100644
index 0000000..062680a
--- /dev/null
+++ b/pysits/docs/content/impute_median.md
@@ -0,0 +1,9 @@
+Remove NA using median
+
+Remove NA using median
+
+Args:
+    data (list[float] | SITSMatrix): A time series vector or matrix.
+
+Returns:
+    R: A set of filtered time series using the imputation function.
diff --git a/pysits/docs/content/plot.md b/pysits/docs/content/plot.md
new file mode 100644
index 0000000..8bdb850
--- /dev/null
+++ b/pysits/docs/content/plot.md
@@ -0,0 +1,112 @@
+Plot sits objects.
+
+Unified plotting function that dispatches on the type of the object passed
+as `x`. It mirrors the many `plot` methods of the R `sits` package,
+covering data cubes (raster, SAR, DEM, vector, RGB), probability and
+uncertainty cubes, variance cubes, classified images, time series patterns
+and predictions, machine learning / deep learning models, clustering and
+self-organizing map (SOM) results, accuracy tables, and t-SNE / embedding
+visualizations. The set of accepted keyword arguments depends on the type
+of object being plotted.
+
+Args:
+    x (SITSCubeModel | SITSTimeSeriesModel | SITSTimeSeriesPatternsModel | SITSMachineLearningMethod | SITSConfusionMatrix): Object to be
+        plotted. Supported objects include classified raster images,
+        classified segments, digital elevation model cubes, multi-year
+        land use/cover embedding predictions, sample distances, class
+        temporal patterns, probability cubes, raster, SAR, and vector
+        data cubes, confusion matrices / accuracy metrics, dendrograms,
+        trained models, time series predictions, t-SNE projections, SOM
+        results, uncertainty cubes, and variance cubes.
+    y: Ignored. Present for compatibility with the generic `plot`.
+    band (str): Band used for plotting a single-band (grey scale) image.
+        Applies to raster, SAR, DEM, and vector cubes, and to SOM maps.
+    red (str): Band assigned to the red channel of an RGB composite
+        (raster, SAR, and vector cubes).
+    green (str): Band assigned to the green channel of an RGB composite.
+    blue (str): Band assigned to the blue channel of an RGB composite.
+    tile (str): Tile to be plotted (data cubes, probability, uncertainty,
+        and variance cubes).
+    dates (list[str]): Dates to be plotted (raster, SAR, and vector
+        cubes).
+    roi (dict | geopandas.GeoDataFrame): Spatial extent (region of
+        interest) to plot, in WGS 84.
+    labels (list[str]): Labels to plot (probability and variance cubes).
+    bands (list[str]): Bands to be viewed (patterns and time series
+        predictions).
+    legend (dict): Associates labels to colors, or a legend specification
+        for SOM plots.
+    legend_position (str): Where to place the legend (typically "inside"
+        or "outside", with defaults varying by plot type).
+    legend_title (str): Title of the legend (probability and variance
+        cubes).
+    palette (str): An RColorBrewer or "cols4all" (or HCL) palette used
+        for color mapping.
+    rev (bool): Whether to reverse the color order in the palette.
+    scale (float): Relative scale of plot text and map (typically 0.4 to
+        1.0).
+    quantile (float): Minimum quantile to plot (probability and variance
+        cubes).
+    first_quantile (float): First quantile for stretching images.
+    last_quantile (float): Last quantile for stretching images.
+    max_cog_size (int): Maximum size of COG (Cloud Optimized GeoTIFF)
+        overviews, in lines/columns or pixels.
+    seg_color (str): Color used to draw segment boundaries (vector cubes).
+    line_width (float): Line width used to draw segment boundaries
+        (vector cubes).
+    type (str): Type of plot; meaning depends on the object. For accuracy
+        objects it is "confusion_matrix" or "metrics"; for variance cubes
+        it is "map" or "hist"; for SOM maps it is "codes" or "mapping".
+    cluster: Cluster object produced by `sits_cluster_dendro`, used when
+        plotting a dendrogram.
+    cutree_height (float): Height at which to draw a dashed horizontal
+        line indicating where the dendrogram is cut.
+    name_cluster (str): Cluster to plot (SOM cluster evaluation).
+    title (str): Title of the plot (SOM cluster evaluation).
+    year_grid (bool): Whether to plot patterns as a grid of panels with
+        labels as columns and years as rows. Defaults to False.
+    tree_idx (int): Index of the tree to be plotted for an XGBoost model.
+    plot_embedding (str): For embedding predictions, either "none" (plot
+        only predicted class intervals) or "area" (overlay a smoothed
+        vertical embedding profile per year).
+    stretch (tuple[float, float]): For embedding plots, lower/upper
+        quantiles used to stretch embedding values before plotting.
+    class_alpha (float): Transparency of class polygons in embedding plots
+        (0-1).
+    area_alpha (float): Transparency of the embedding area in embedding
+        plots (0-1).
+    area_width (float): Horizontal width fraction of the embedding area.
+    area_spar (float): Smoothing parameter for the embedding area spline.
+    **kwargs (dict): Further specifications passed to the underlying plot.
+
+Returns:
+    None: A plot is produced. Depending on the input type this may be a
+    color map of classified pixels, an RGB or grey-scale image, a
+    probability or uncertainty map, a variance map (optionally with
+    segment overlays), a dendrogram, a confusion matrix, a SOM map, a
+    model diagnostic plot, or a plot for patterns, predictions,
+    embeddings, and t-SNE projections. Some methods are called only for
+    their side effect of drawing the plot.
+
+Notes:
+    The `roi` argument can be defined as a `dict` giving the spatial
+    extent (for example with `lon_min`, `lon_max`, `lat_min`, `lat_max`),
+    a `geopandas.GeoDataFrame`, or another spatial specification accepted
+    by `sits`. Vector cube plots overlay the segments produced by
+    `sits_segment` on top of the raster image; their appearance is
+    controlled by `seg_color` and `line_width`.
+
+Examples:
+    from pysits import *
+
+    # Plot a set of time series patterns
+    patterns = sits_patterns(cerrado_2classes)
+    plot(patterns)
+
+    # Train a random forest model and plot variable importance
+    rfor_model = sits_train(samples_modis_ndvi, ml_method=sits_rfor())
+    plot(rfor_model)
+
+    # Plot a SOM map produced from a set of samples
+    som_map = sits_som_map(samples_modis_ndvi)
+    plot(som_map)
diff --git a/pysits/docs/content/sits_accuracy.md b/pysits/docs/content/sits_accuracy.md
new file mode 100644
index 0000000..38dc63e
--- /dev/null
+++ b/pysits/docs/content/sits_accuracy.md
@@ -0,0 +1,80 @@
+Assess classification accuracy
+
+This function calculates the accuracy of the classification result. The input
+is either a set of classified time series or a classified data cube. Classified
+time series are produced by `sits_classify`. Classified images are generated
+using `sits_classify` followed by `sits_label_classification`.
+For a set of time series, `sits_accuracy` creates a confusion matrix and
+calculates the resulting statistics using package `caret`. For a classified
+image, the function uses an area-weighted technique proposed by Olofsson et al.
+according to references [1-3] to produce reliable accuracy estimates at 95%
+confidence level. In both cases, it provides an accuracy assessment of the
+classified, including Overall Accuracy, Kappa, User's Accuracy, Producer's
+Accuracy and error matrix (confusion matrix).
+
+Args:
+    data (SITSCubeModel | SITSTimeSeriesModel): Either a data cube with
+        classified images or a set of time series.
+    prediction_attr (str): Name of the column of the segments object that
+        contains the predicted values (only for vector class cubes).
+    reference_attr (str): Name of the column of the segments object that
+        contains the reference values (only for vector class cubes).
+    validation (str | pathlib.Path | pandas.DataFrame | geopandas.GeoDataFrame | SITSTimeSeriesModel):
+        Samples for validation (see below). Only required when data is a
+        raster class cube.
+    method (str): Either 'olofsson' or 'pixel' to compute accuracy (only
+        for raster class cubes).
+    **kwargs (dict): Specific parameters.
+
+Returns:
+    SITSData: The error_matrix, the class_areas, the unbiased estimated
+    areas, the standard error areas, confidence interval 95 and the accuracy
+    (user, producer, and overall), or `None` if the data is empty. The result
+    can be visualized directly on the screen.
+
+Notes:
+    The `validation` data needs to contain the following columns: "latitude",
+    "longitude", "start_date", "end_date", and "label". It can be either a path
+    to a CSV file, a `SITSTimeSeriesModel`, a `pandas.DataFrame`, or a
+    `geopandas.GeoDataFrame`.
+    When `validation` is a `geopandas.GeoDataFrame`, the columns "latitude"
+    and "longitude" are not required as the locations are extracted from the
+    geometry column. The `centroid` is calculated before extracting the
+    location values for any geometry type.
+
+Examples:
+    from pysits import *
+
+    # show accuracy for a set of samples
+    train_data = sits_sample(samples_modis_ndvi, frac=0.5)
+    test_data = sits_sample(samples_modis_ndvi, frac=0.5)
+    rfor_model = sits_train(train_data, sits_rfor())
+    points_class = sits_classify(
+        data=test_data, ml_model=rfor_model
+    )
+    acc = sits_accuracy(points_class)
+
+    # show accuracy for a data cube classification
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    import tempfile
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # label the probability cube
+    label_cube = sits_label_classification(
+        probs_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # obtain the ground truth for accuracy assessment
+    ground_truth = r_package_dir("extdata/samples/samples_sinop_crop.csv", package="sits")
+    # make accuracy assessment
+    as_ = sits_accuracy(label_cube, validation=ground_truth)
diff --git a/pysits/docs/content/sits_accuracy_summary.md b/pysits/docs/content/sits_accuracy_summary.md
new file mode 100644
index 0000000..877c47e
--- /dev/null
+++ b/pysits/docs/content/sits_accuracy_summary.md
@@ -0,0 +1,11 @@
+Print accuracy summary
+
+Adaptation of the caret::print.confusionMatrix method for the more common
+usage in Earth Observation.
+
+Args:
+    x (SITSConfusionMatrix): accuracy object to summarize.
+    digits (int): number of significant digits when printed.
+
+Returns:
+    SITSData: called for side effects.
diff --git a/pysits/docs/content/sits_add_base_cube.md b/pysits/docs/content/sits_add_base_cube.md
new file mode 100644
index 0000000..a3b3bd2
--- /dev/null
+++ b/pysits/docs/content/sits_add_base_cube.md
@@ -0,0 +1,55 @@
+Add base maps to a time series data cube
+
+This function add base maps to time series data cube. Base maps have
+information that is stable in time (e.g, DEM) which provide relevant
+information for modelling and classification.
+To add a base cube to an existing data cube, they should share the same
+sensor, resolution, bounding box, timeline, and have different bands.
+
+Args:
+    cube1 (SITSCubeModel): Data cube.
+    cube2 (SITSCubeModel): Data cube with base information.
+
+Returns:
+    SITSCubeModel: a merged data cube with the inclusion of base
+        information.
+
+Examples:
+    from pysits import *
+    import tempfile
+    import os
+
+    s2_cube = sits_cube(
+        source="MPC",
+        collection="SENTINEL-2-L2A",
+        tiles="18HYE",
+        bands=["B8A", "CLOUD"],
+        start_date="2022-01-01",
+        end_date="2022-03-31"
+    )
+    output_dir = os.path.join(tempfile.gettempdir(), "reg")
+    if not os.path.exists(output_dir):
+        os.makedirs(output_dir)
+    dem_cube = sits_cube(
+        source="MPC",
+        collection="COP-DEM-GLO-30",
+        tiles="18HYE",
+        bands="ELEVATION"
+    )
+    s2_reg = sits_regularize(
+        cube=s2_cube,
+        period="P1M",
+        res=240,
+        output_dir=output_dir,
+        multicores=2,
+        memsize=4
+    )
+    dem_reg = sits_regularize(
+        cube=dem_cube,
+        res=240,
+        tiles="18HYE",
+        output_dir=output_dir,
+        multicores=2,
+        memsize=4
+    )
+    s2_reg = sits_add_base_cube(s2_reg, dem_reg)
diff --git a/pysits/docs/content/sits_apply.md b/pysits/docs/content/sits_apply.md
new file mode 100644
index 0000000..473066d
--- /dev/null
+++ b/pysits/docs/content/sits_apply.md
@@ -0,0 +1,104 @@
+Apply a function on a set of time series
+
+Apply a named expression to a set of time series or a data cube to be
+evaluated and generate new bands (indices). In the case of data cubes,
+it creates a new band in `output_dir`.
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): valid time series or data
+        cube.
+    window_size (int): an odd number representing the size of the sliding
+        window of kernel functions used in expressions (for a list of
+        supported kernel functions, please see details).
+    memsize (int): memory available for classification (in GB).
+    multicores (int): number of cores to be used for classification.
+    normalized (bool): does the expression produce a normalized band?
+    output_dir (str | pathlib.Path): directory where files will be saved.
+    progress (bool): show progress bar?
+    **kwargs (dict): named expressions to be evaluated (see details).
+
+Returns:
+    SITSFrame: time series or data cube with new bands, produced according
+        to the requested expression.
+
+Notes:
+    The main `sits` classification workflow has the following steps:
+    1. `sits_cube`: selects a ARD image collection from a cloud provider.
+    2. `sits_cube_copy`: copies an ARD image collection from a cloud provider
+       to a local directory for faster processing.
+    3. `sits_regularize`: create a regular data cube from an ARD image
+       collection.
+    4. `sits_apply`: create new indices by combining bands of a regular data
+       cube (optional).
+    5. `sits_get_data`: extract time series from a regular data cube based on
+       user-provided labelled samples.
+    6. `sits_train`: train a machine learning model based on image time series.
+    7. `sits_classify`: classify a data cube using a machine learning model and
+       obtain a probability cube.
+    8. `sits_smooth`: post-process a probability cube using a spatial smoother
+       to remove outliers and increase spatial consistency.
+    9. `sits_label_classification`: produce a classified map by selecting the
+       label with the highest probability from a smoothed cube.
+    `sits_apply()` allows any valid R expression to compute new bands. Use R
+    syntax to pass an expression to this function. Besides arithmetic
+    operators, you can use virtually any R function that can be applied to
+    elements of a matrix (functions that are unaware of matrix sizes, e.g.
+    `sqrt()`, `sin()`, `log()`).
+    Examples of valid expressions:
+    1. `NDVI = (B08 - B04) / (B08 + B04)` for Sentinel-2 images.
+    2. `EVI = 2.5 * (B05 \04) / (B05 + 6 * B04 \7.5 * B02 + 1)` for
+       Landsat-8/9 images.
+    3. `VV_VH_RATIO = VH/VV` for Sentinel-1 images. In this case, set the
+       `normalized` parameter to `False`.
+    4. `VV_DB = 10 * log10(VV)` to convert Sentinel-1 RTC images available in
+       Planetary Computer to decibels. Also, set the `normalized` parameter to
+       `False`.
+    `sits_apply()` also accepts a predefined set of kernel functions (see
+    below) that can be applied to pixels considering its neighborhood. The
+    function considers a neighborhood of a pixel as a set of pixels equidistant
+    to it (including itself). This neighborhood forms a square window (also
+    known as kernel) around the central pixel (Moore neighborhood). Users can
+    set the `window_size` parameter to adjust the size of the kernel window.
+    The image is conceptually mirrored at the edges so that neighborhood
+    including a pixel outside the image is equivalent to take the 'mirrored'
+    pixel inside the edge.
+    `sits_apply()` applies a function to the kernel and its result is assigned
+    to a corresponding central pixel on a new matrix. The kernel slides
+    throughout the input image and this process generates an entire new matrix,
+    which is returned as a new band to the cube. The kernel functions ignores
+    any `NA` values inside the kernel window. If all pixels in the window are
+    `NA` the result will be `NA`.
+    By default, the indexes generated by `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 `normalized` parameter is
+    `False`, no scaling factor will be applied and the index will be saved as
+    FLT4S (signed float) and the values will vary between -3.4e+38 and
+    3.4e+38.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # get a time series
+    # Apply a normalization function
+    point2 = sits_apply(
+        sits_select(point_mt_6bands, "NDVI"),
+        NDVI_norm="(NDVI - min(NDVI)) / (max(NDVI) - min(NDVI))"
+    )
+
+    # Example of generation texture band with variance
+    # Create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # Generate a texture images with variance in NDVI images
+    cube_texture = sits_apply(
+        data=cube,
+        NDVITEXTURE="w_median(NDVI)",
+        window_size=5,
+        output_dir=tempfile.mkdtemp()
+    )
diff --git a/pysits/docs/content/sits_as_geopandas.md b/pysits/docs/content/sits_as_geopandas.md
new file mode 100644
index 0000000..15466e7
--- /dev/null
+++ b/pysits/docs/content/sits_as_geopandas.md
@@ -0,0 +1,28 @@
+Return time series or a data cube as a `geopandas.GeoDataFrame`.
+
+Converts time series or a data cube to a `geopandas.GeoDataFrame`.
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): time series or data
+            cube.
+    crs (str): input coordinate reference system.
+    as_crs (str): output coordinate reference system.
+    **kwargs (dict): additional parameters.
+
+Returns:
+    SITSFrame: point or polygon geometry.
+
+Examples:
+    from pysits import *
+
+    # convert sits tibble to a geopandas object (point)
+    geo_object = sits_as_geopandas(cerrado_2classes)
+
+    # convert sits cube to a geopandas object (polygon)
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    geo_object = sits_as_geopandas(cube)
diff --git a/pysits/docs/content/sits_bands.md b/pysits/docs/content/sits_bands.md
new file mode 100644
index 0000000..140779d
--- /dev/null
+++ b/pysits/docs/content/sits_bands.md
@@ -0,0 +1,30 @@
+Get the names of the bands
+
+Finds the names of the bands of a set of time series or of a data cube
+
+Args:
+    x (SITSTimeSeriesModel | SITSCubeModel): time series or data cube.
+    value (list[str]): new value for the bands.
+
+Returns:
+    list: the names of the bands.
+
+Examples:
+    from pysits import *
+
+    # Create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # Get the bands from a data cube
+    bands = sits_bands(cube)
+    # Get the bands from a sits tibble
+    bands = sits_bands(samples_modis_ndvi)
+    # Get the bands from patterns
+    bands = sits_bands(sits_patterns(samples_modis_ndvi))
+    # Get the bands from ML model
+    rf_model = sits_train(samples_modis_ndvi, sits_rfor())
+    bands = sits_bands(rf_model)
diff --git a/pysits/docs/content/sits_barlow_twins.md b/pysits/docs/content/sits_barlow_twins.md
new file mode 100644
index 0000000..dccb95a
--- /dev/null
+++ b/pysits/docs/content/sits_barlow_twins.md
@@ -0,0 +1,65 @@
+Barlow Twins encoder for image time series
+
+Supervised pre-training using the Barlow Twins loss and a torch encoder. Two
+time series with the same class label are passed through a shared encoder +
+projector. The Barlow Twins loss makes the cross-correlation matrix of the two
+views' embeddings close to the identity: the diagonal -> 1 (invariance) and the
+off-diagonal -> 0 (redundancy reduction). No negatives are required.
+The function can be used in two ways:
+- If `samples` is provided, it trains immediately and returns an encoder-ready
+  model object (see Value).
+- If `samples = None`, it returns a training function with signature
+  `function(samples)` that can be passed to `sits_pre_train` or called later.
+
+Args:
+    samples (SITSTimeSeriesModel): A set of sample time series. If `None`
+        (default), returns a training function. If provided, triggers
+        immediate training. Base data samples (e.g., `sits_base`) are not
+        supported.
+    embedding_dim (int): Dimensionality of the encoder embedding (exported
+        features). Default: 64.
+    proj_dim (int): Dimensionality of the projector head used only during
+        pre-training. Default: 256.
+    bt_lambda (float): Weight of the redundancy-reduction (off-diagonal)
+        term in the Barlow Twins loss. Default: 5e-3.
+    num_pairs (int | None): Total number of pairs to form per epoch. When
+        `None` (default), one pair is formed for every sample in the
+        training split.
+    encoder_model (SITSMachineLearningMethod): Deep learning method that
+        takes time series as input and produces latent representations that
+        are used to compute the loss function (suggested options:
+        `sits_tempcnn()`, `sits_lighttae()`, `sits_resnet()`). Default:
+        `sits_tempcnn()`.
+    epochs (int): Maximum number of training epochs.
+    batch_size (int): Batch size for training. Larger values improve the
+        Barlow Twins cross-correlation estimate. Default: 128.
+    validation_split (float): Fraction of samples held out for validation
+        loss monitoring, in the range (0, 1).
+    optimizer: A `torch` optimizer constructor (default:
+        `torch::optim_adamw`).
+    opt_hparams (dict): Optimizer hyperparameters. Common entries: `lr`,
+        `eps`, `weight_decay`.
+    lr_decay_epochs (int): Step size (in epochs) for LR decay.
+    lr_decay_rate (float): Multiplicative LR decay factor.
+    patience (int): Early-stopping patience (epochs without improvement).
+    min_delta (float): Minimum improvement required to reset the patience
+        counter.
+    verbose (bool): Print training progress?
+    seed (int): Random seed for reproducibility.
+
+Returns:
+    R: If `samples = None`, a training function with signature
+    `function(samples)` that trains a Barlow Twins model and returns a
+    pretrained encoder. If `samples` is provided, the result of applying the
+    training function to `samples` directly.
+
+Examples:
+    from pysits import *
+
+    model = sits_pre_train(
+        samples_modis_ndvi,
+        sits_barlow_twins(
+            embedding_dim=32,
+            epochs=20
+        )
+    )
diff --git a/pysits/docs/content/sits_bbox.md b/pysits/docs/content/sits_bbox.md
new file mode 100644
index 0000000..9f24243
--- /dev/null
+++ b/pysits/docs/content/sits_bbox.md
@@ -0,0 +1,33 @@
+Get the bounding box of the data
+
+Obtain a vector of limits (either on lat/long for time series or in
+projection coordinates in the case of cubes)
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): samples or data cube.
+    crs (str): CRS of the time series.
+    as_crs (str): CRS to project the resulting bounding box.
+    **kwargs (dict): parameters for specific types.
+
+Returns:
+    SITSFrame: the bounding box.
+
+Notes:
+    Time series in `sits` are associated with lat/long values in WGS84,
+    while each data cube is associated to a cartographic projection. To
+    obtain the bounding box of a data cube in a different projection than
+    the original, use the `as_crs` parameter.
+
+Examples:
+    from pysits import *
+
+    # get the bbox of a set of samples
+    sits_bbox(samples_modis_ndvi)
+    # get the bbox of a cube in WGS84
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    sits_bbox(cube, as_crs="EPSG:4326")
diff --git a/pysits/docs/content/sits_classify.md b/pysits/docs/content/sits_classify.md
new file mode 100644
index 0000000..566cd8f
--- /dev/null
+++ b/pysits/docs/content/sits_classify.md
@@ -0,0 +1,107 @@
+Classify a set of time series or a data cube.
+
+This function applies a machine learning model (trained by `sits_train`)
+to classify time series or data cubes. Its behavior depends on the type
+of the input `data`:
+
+- Set of time series (`SITSTimeSeriesModel`): the output is the same set
+  of time series with an additional `predicted` column containing the
+  assigned labels for each point.
+- Regular raster cube (`SITSCubeModel`): the output is a probability cube
+  with the same tiles as the input. Each tile is a multiband image where
+  each band contains the probability that a pixel belongs to a given
+  class.
+- Segmented (vector) data cube (produced by `sits_segment`): the temporal
+  model is applied to produce pixel-level probabilities, and the
+  associated vector support (`vector_info`) is preserved in the output.
+  The result is a probability cube with vector support, which can then be
+  passed to `sits_label_classification` for segment-based labeling.
+  Segment-level aggregation is no longer performed by `sits_classify`;
+  use `sits_label_classification` to aggregate pixel probabilities inside
+  segments and assign classes.
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): input to classify. Either
+        a set of time series, a regular raster data cube, or a segmented
+        vector data cube.
+    ml_model (SITSMachineLearningMethod): model trained by `sits_train`.
+    roi (dict | str | pathlib.Path | geopandas.GeoDataFrame): region of
+        interest, either a `geopandas.GeoDataFrame`, a shapefile, or a
+        `dict` in WGS 84 with named XY values (`xmin`, `xmax`, `ymin`,
+        `ymax`) or named lat/long values (`lon_min`, `lat_min`,
+        `lon_max`, `lat_max`). Applies to raster and vector cubes.
+    exclusion_mask (geopandas.GeoDataFrame | str | pathlib.Path): areas
+        to be excluded from the classification process, defined by a
+        `geopandas.GeoDataFrame` or a shapefile. Applies to raster and
+        vector cubes.
+    impute_fn: imputation function to remove NA.
+    start_date (str): starting date for the classification (in YYYY-MM-DD
+        format). Applies to raster and vector cubes.
+    end_date (str): ending date for the classification (in YYYY-MM-DD
+        format). Applies to raster and vector cubes.
+    memsize (int): memory available for classification in GB (min = 1,
+        max = 16384). Applies to raster and vector cubes.
+    multicores (int): number of cores to be used for classification
+        (min = 1, max = 2048).
+    gpu_memory (int): memory available in GPU in GB (default = 4).
+    batch_size (int): batch size for GPU classification.
+    block_size (dict): size of the block read and written by each worker,
+        with `[nrows, ncols]`. Default is `None`, which computes an
+        optimal block size from `memsize`, `multicores` and the internal
+        block size of the raster files. Applies to raster and vector
+        cubes.
+    output_dir (str | pathlib.Path): directory for output file. Applies
+        to raster and vector cubes.
+    version (str): version of the output. Applies to raster and vector
+        cubes.
+    n_sam_pol (int): deprecated. Segment-level classification is no longer
+        performed by `sits_classify`. Use `sits_label_classification` for
+        segment-based labeling. Applies to vector cubes.
+    verbose (bool): print information about processing time? Applies to
+        raster and vector cubes.
+    progress (bool): show progress bar?
+    **kwargs (dict): other parameters for specific functions.
+
+Returns:
+    SITSCubeModel: for a set of time series, a `SITSTimeSeriesModel` with
+        predicted labels for each point. For a regular raster cube, a data
+        cube with probabilities for each class. For a segmented vector
+        cube, a probability cube with associated vector support that
+        contains pixel-level probabilities and preserves `vector_info` for
+        segment-based labeling.
+
+Examples:
+    from pysits import *
+
+    # Example 1: classify a set of time series
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+
+    # classify a point
+    point_ndvi = sits_select(point_mt_6bands, bands=["NDVI"])
+    point_class = sits_classify(data=point_ndvi, ml_model=rfor_model)
+    plot(point_class)
+
+    # Example 2: classify a raster cube
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube,
+        ml_model=rfor_model,
+        output_dir="./tempdir"
+    )
+    plot(probs_cube)
+
+    # label the probability cube
+    label_cube = sits_label_classification(
+        probs_cube,
+        output_dir="./tempdir"
+    )
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_clean.md b/pysits/docs/content/sits_clean.md
new file mode 100644
index 0000000..bec6122
--- /dev/null
+++ b/pysits/docs/content/sits_clean.md
@@ -0,0 +1,58 @@
+Cleans a classified map using a local window
+
+Applies a modal function to clean up possible noisy pixels keeping the most
+frequently values within the neighborhood. In a tie, the first value of the
+vector is considered. Modal functions applied to classified cubes are useful to
+remove salt-and-pepper noise in the result.
+
+Args:
+    cube (SITSCubeModel): Classified data cube.
+    window_size (int): An odd integer representing the size of the sliding
+        window of the modal function (min = 1, max = 15).
+    memsize (int): Memory available for classification in GB (min = 1,
+        max = 16384).
+    multicores (int): Number of cores to be used for classification (min = 1,
+        max = 2048).
+    output_dir (str | pathlib.Path): Valid directory for output file.
+    version (str): Version of the output file.
+    progress (bool): Show progress bar?
+    **kwargs (dict): Specific parameters for specialised functions.
+
+Returns:
+    SITSCubeModel: A classified map.
+
+Notes:
+    The `sits_clean` function is useful to further remove classification noise
+    which has not been detected by `sits_smooth`. It improves the spatial
+    consistency of the classified maps.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    rf_model = sits_train(samples_modis_ndvi, ml_method=sits_rfor)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube,
+        ml_model=rf_model,
+        output_dir=tempfile.gettempdir()
+    )
+    # label the probability cube
+    label_cube = sits_label_classification(
+        probs_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # apply a mode function in the labelled cube
+    clean_cube = sits_clean(
+        cube=label_cube,
+        window_size=5,
+        output_dir=tempfile.gettempdir(),
+        multicores=1
+    )
diff --git a/pysits/docs/content/sits_cluster_clean.md b/pysits/docs/content/sits_cluster_clean.md
new file mode 100644
index 0000000..2c327b3
--- /dev/null
+++ b/pysits/docs/content/sits_cluster_clean.md
@@ -0,0 +1,22 @@
+Removes labels that are minority in each cluster.
+
+Takes a set of time series that has additional `cluster` information
+produced by `sits_cluster_dendro()` and removes labels that are minority
+in each cluster.
+
+Args:
+    samples (SITSTimeSeriesModel): Set of time series with additional
+        cluster information produced by `sits_cluster_dendro()`.
+
+Returns:
+    SITSTimeSeriesModel: Set of time series.
+
+Examples:
+    from pysits import *
+
+    clusters = sits_cluster_dendro(cerrado_2classes)
+    freq1 = sits_cluster_frequency(clusters)
+    print(freq1)
+    clean_clusters = sits_cluster_clean(clusters)
+    freq2 = sits_cluster_frequency(clean_clusters)
+    print(freq2)
diff --git a/pysits/docs/content/sits_cluster_dendro.md b/pysits/docs/content/sits_cluster_dendro.md
new file mode 100644
index 0000000..db7e4fe
--- /dev/null
+++ b/pysits/docs/content/sits_cluster_dendro.md
@@ -0,0 +1,47 @@
+Find clusters in time series samples
+
+These functions support hierarchical agglomerative clustering in sits. They
+provide support from creating a dendrogram and using it for cleaning samples.
+`sits_cluster_dendro()` takes time series and produces a `SITSTimeSeriesModel`
+with an added "cluster" column. The function first calculates a dendrogram and
+obtains a validity index for best clustering using the adjusted Rand Index.
+After cutting the dendrogram using the chosen validity index, it assigns a
+cluster to each sample.
+`sits_cluster_frequency()` computes the contingency table between labels and
+clusters and produces a matrix. Its input is produced by
+`sits_cluster_dendro()`.
+`sits_cluster_clean()` takes time series that have an additional `cluster`
+produced by `sits_cluster_dendro()` and removes labels that are minority in
+each cluster.
+
+Args:
+    samples (SITSTimeSeriesModel): input set of time series.
+    bands (list[str]): bands to be used in the clustering.
+    dist_method (str): one of the supported distances. "dtw": DTW with a
+        Sakoe-Chiba constraint. "dtw2": DTW with L2 norm and Sakoe-Chiba
+        constraint. "dtw_basic": A faster DTW with less functionality. "lbk":
+        Keogh's lower bound for DTW. "lbi": Lemire's lower bound for DTW.
+    linkage (str): agglomeration method to be used. One of "ward.D",
+        "ward.D2", "single", "complete", "average", "mcquitty", "median" or
+        "centroid".
+    k (int): desired number of clusters (overrides default value).
+    palette (str): color palette as per `grDevices::hcl.pals()` function.
+    **kwargs (dict): additional parameters to be passed to
+        dtwclust::tsclust() function.
+
+Returns:
+    SITSTimeSeriesModel: time series with an added "cluster" column.
+
+Notes:
+    Please refer to the sits documentation available in
+    https://e-sensing.github.io/sitsbook/ for detailed examples.
+
+Examples:
+    from pysits import *
+
+    # default
+    clusters = sits_cluster_dendro(cerrado_2classes)
+    # with parameters
+    clusters = sits_cluster_dendro(cerrado_2classes,
+        bands="NDVI", k=5
+    )
diff --git a/pysits/docs/content/sits_cluster_frequency.md b/pysits/docs/content/sits_cluster_frequency.md
new file mode 100644
index 0000000..68a17fc
--- /dev/null
+++ b/pysits/docs/content/sits_cluster_frequency.md
@@ -0,0 +1,18 @@
+Show label frequency in each cluster produced by dendrogram analysis
+
+Show label frequency in each cluster produced by dendrogram analysis
+
+Args:
+    samples (SITSTimeSeriesModel): input set of time series with
+        additional cluster information produced by
+        `sits_cluster_dendro()`.
+
+Returns:
+    SITSTable: frequencies of labels in clusters.
+
+Examples:
+    from pysits import *
+
+    clusters = sits_cluster_dendro(cerrado_2classes)
+    freq = sits_cluster_frequency(clusters)
+    print(freq)
diff --git a/pysits/docs/content/sits_colors.md b/pysits/docs/content/sits_colors.md
new file mode 100644
index 0000000..8f1dae8
--- /dev/null
+++ b/pysits/docs/content/sits_colors.md
@@ -0,0 +1,22 @@
+Function to retrieve sits color table
+
+Returns the default color table.
+
+Args:
+    legend (str): One of the accepted legends in sits.
+
+Returns:
+    SITSFrame: color names and values.
+
+Notes:
+    SITS has a predefined color palette with 238 class names. These colors are
+    grouped by typical legends used by the Earth observation community, which
+    include "IGBP", "UMD", "ESA_CCI_LC", and "WORLDCOVER". Use
+    `sits_colors_show` to see a specific palette. The default color table can
+    be extended using `sits_colors_set`.
+
+Examples:
+    from pysits import *
+
+    # return the names of all colors supported by SITS
+    sits_colors()
diff --git a/pysits/docs/content/sits_colors_qgis.md b/pysits/docs/content/sits_colors_qgis.md
new file mode 100644
index 0000000..68d1132
--- /dev/null
+++ b/pysits/docs/content/sits_colors_qgis.md
@@ -0,0 +1,34 @@
+Function to save color table as QML style for data cube
+
+Saves a color table associated to a classified data cube as a QGIS
+style file
+
+Args:
+    cube (SITSCubeModel): a classified data cube.
+    file (str | pathlib.Path): a QGIS style file to be written to.
+
+Returns:
+    None: called for side effects.
+
+Examples:
+    from pysits import *
+
+    data_dir = r_package_dir("extdata/raster/classif", package="sits")
+    ro_class = sits_cube(
+        source="MPC",
+        collection="SENTINEL-2-L2A",
+        data_dir=data_dir,
+        parse_info=[
+            "X1", "X2", "tile", "start_date", "end_date",
+            "band", "version"
+        ],
+        bands="class",
+        labels={
+            "1": "Clear_Cut_Burned_Area",
+            "2": "Clear_Cut_Bare_Soil",
+            "3": "Clear_Cut_Vegetation",
+            "4": "Forest"
+        }
+    )
+    qml_file = "/tmp/qgis.qml"
+    sits_colors_qgis(ro_class, qml_file)
diff --git a/pysits/docs/content/sits_colors_reset.md b/pysits/docs/content/sits_colors_reset.md
new file mode 100644
index 0000000..9ca281b
--- /dev/null
+++ b/pysits/docs/content/sits_colors_reset.md
@@ -0,0 +1,12 @@
+Function to reset sits color table
+
+Resets the color table
+
+Returns:
+    None: called for side effects.
+
+Examples:
+    from pysits import *
+
+    # reset the default colors supported by SITS
+    sits_colors_reset()
diff --git a/pysits/docs/content/sits_colors_set.md b/pysits/docs/content/sits_colors_set.md
new file mode 100644
index 0000000..82cd65c
--- /dev/null
+++ b/pysits/docs/content/sits_colors_set.md
@@ -0,0 +1,69 @@
+Function to set sits color table
+
+Includes new colors in the SITS color sets. If the colors exist, replace them
+with the new HEX value. Optionally, the new colors can be associated to a
+legend. In this case, the new legend name should be informed. The colors
+parameter should be a `pandas.DataFrame` with name and HEX code. Colour
+names should be one character string only. Composite names need to be combined
+with underscores (e.g., use "Snow_and_Ice" and not "Snow and Ice").
+This function changes the global sits color table and the global set of sits
+color legends. To undo these effects, please use "sits_colors_reset()".
+
+Args:
+    colors (pandas.DataFrame): New color table with name and HEX code.
+    legend (str): Legend associated to the color table (optional).
+
+Returns:
+    SITSFrame: A modified sits color table (invisible).
+
+Examples:
+    from pysits import *
+    import pandas as pd
+
+    # Define a color table based on the Anderson Land Classification System
+    us_nlcd = pd.DataFrame({
+        "name": [
+            "Urban_Built_Up",
+            "Agricultural_Land",
+            "Rangeland",
+            "Forest_Land",
+            "Water",
+            "Wetland",
+            "Barren_Land",
+            "Tundra",
+            "Snow_and_Ice",
+        ],
+        "color": [
+            "#85929E",
+            "#F0B27A",
+            "#F1C40F",
+            "#27AE60",
+            "#2980B9",
+            "#D4E6F1",
+            "#FDEBD0",
+            "#EBDEF0",
+            "#F7F9F9",
+        ],
+    })
+
+    # Load the color table into `sits`
+    sits_colors_set(colors=us_nlcd, legend="US_NLCD")
+
+    # Show the new color table used by sits
+    sits_colors_show("US_NLCD")
+
+    # Change colors in the sits global color table
+    # First show the default colors for the UMD legend
+    sits_colors_show("UMD")
+    # Then change some colors associated to the UMD legend
+    mycolors = pd.DataFrame({
+        "name": ["Savannas", "Grasslands"],
+        "color": ["#F8C471", "#ABEBC6"],
+    })
+    sits_colors_set(colors=mycolors)
+    # Notice that the UMD colors change
+    sits_colors_show("UMD")
+    # Reset the color table
+    sits_colors_reset()
+    # Show the default colors for the UMD legend
+    sits_colors_show("UMD")
diff --git a/pysits/docs/content/sits_colors_show.md b/pysits/docs/content/sits_colors_show.md
new file mode 100644
index 0000000..7e0a298
--- /dev/null
+++ b/pysits/docs/content/sits_colors_show.md
@@ -0,0 +1,14 @@
+Function to show colors in SITS
+
+Shows the default SITS colors
+
+Args:
+    legend (str): One of the accepted legends in sits.
+    font_family (str): A font family loaded in SITS.
+
+Returns:
+    None: called for side effects.
+
+Examples:
+    # show the colors supported by SITS
+    sits_colors_show()
diff --git a/pysits/docs/content/sits_combine_predictions.md b/pysits/docs/content/sits_combine_predictions.md
new file mode 100644
index 0000000..308acdf
--- /dev/null
+++ b/pysits/docs/content/sits_combine_predictions.md
@@ -0,0 +1,72 @@
+Estimate ensemble prediction based on list of probs cubes
+
+Calculate an ensemble predictor based a list of probability cubes. The function
+combines the output of two or more models to derive a weighted average. The
+supported types of ensemble predictors are 'average' and 'uncertainty'. In the
+latter case, the uncertainty cubes need to be provided using param
+`uncert_cubes`.
+
+Args:
+    cubes (list[SITSCubeModel]): List of probability data cubes.
+    type (str): Method to measure uncertainty. One of "average" or
+        "uncertainty".
+    weights (list[float]): Weights for averaging.
+    memsize (int): Memory available for classification in GB (min = 1,
+        max = 16384).
+    multicores (int): Number of cores to be used for classification
+        (min = 1, max = 2048).
+    output_dir (str | pathlib.Path): Valid directory for output file.
+    version (str): Version of the output.
+    progress (bool): Set progress bar?
+    uncert_cubes (list[SITSCubeModel]): Uncertainty cubes to be used as
+        local weights when type = "uncertainty" is selected.
+    **kwargs (dict): Parameters for specific functions.
+
+Returns:
+    SITSCubeModel: A combined probability cube.
+
+Notes:
+    The distribution of class probabilities produced by machine learning models
+    such as random forest is quite different from that produced by deep
+    learning models such as temporal CNN. Combining the result of two different
+    models is recommended to remove possible bias induced by a single model.
+    By default, the function takes the average of the class probabilities of
+    two or more model results. If desired, users can use the uncertainty
+    estimates for each results to compute the weights for each pixel. In this
+    case, the uncertainties produced by the models for each pixel are used to
+    compute the weights for producing the combined result.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # classify a data cube using rfor model
+    probs_rfor_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir(),
+        version="rfor"
+    )
+    # create an SVM model
+    svm_model = sits_train(samples_modis_ndvi, sits_svm())
+    # classify a data cube using SVM model
+    probs_svm_cube = sits_classify(
+        data=cube, ml_model=svm_model, output_dir=tempfile.gettempdir(),
+        version="svm"
+    )
+    # create a list of predictions to be combined
+    pred_cubes = [probs_rfor_cube, probs_svm_cube]
+    # combine predictions
+    comb_probs_cube = sits_combine_predictions(
+        pred_cubes,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the resulting combined prediction cube
+    plot(comb_probs_cube)
diff --git a/pysits/docs/content/sits_confidence_sampling.md b/pysits/docs/content/sits_confidence_sampling.md
new file mode 100644
index 0000000..e8c73c6
--- /dev/null
+++ b/pysits/docs/content/sits_confidence_sampling.md
@@ -0,0 +1,54 @@
+Suggest high confidence samples to increase the training set.
+
+Suggest points for increasing the training set. These points are labelled with
+high confidence so they can be added to the training set. They need to have a
+satisfactory margin of confidence to be selected. The input is a probability
+cube. For each label, the algorithm finds out location where the machine
+learning model has high confidence in choosing this label compared to all
+others. The algorithm also considers a minimum distance between new labels, to
+minimize spatial autocorrelation effects. This function is best used in the
+following context:
+1. Select an initial set of samples.
+2. Train a machine learning model.
+3. Build a data cube and classify it using the model.
+4. Run a Bayesian smoothing in the resulting probability cube.
+5. Perform confidence sampling.
+The Bayesian smoothing procedure will reduce the classification outliers and
+thus increase the likelihood that the resulting pixels with provide good
+quality samples for each class.
+
+Args:
+    probs_cube (SITSCubeModel): A smoothed probability cube. See
+        `sits_classify` and `sits_smooth`.
+    n (int): Number of suggested points per class.
+    min_margin (float): Minimum margin of confidence to select a sample.
+    sampling_window (int): Window size for collecting points (in pixels).
+        The minimum window size is 10.
+    multicores (int): Number of workers for parallel processing (min = 1,
+        max = 2048).
+    memsize (int): Maximum overall memory (in GB) to run the function.
+    progress (bool): Show progress bar?
+
+Returns:
+    SITSFrame: longitude and latitude in WGS84 with locations which have
+        high uncertainty and meet the minimum distance criteria.
+
+Examples:
+    from pysits import *
+
+    # create a data cube
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # build a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, ml_method=sits_rfor())
+    # classify the cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempdir()
+    )
+    # obtain a new set of samples for active learning
+    # the samples are located in uncertain places
+    new_samples = sits_confidence_sampling(probs_cube)
diff --git a/pysits/docs/content/sits_config.md b/pysits/docs/content/sits_config.md
new file mode 100644
index 0000000..0ddcd82
--- /dev/null
+++ b/pysits/docs/content/sits_config.md
@@ -0,0 +1,28 @@
+Configure parameters for sits package
+
+These functions load and show sits configurations.
+The `sits` package uses a configuration file that contains information
+on parameters required by different functions. This includes
+information about the image collections handled by `sits`.
+`sits_config()` loads the default configuration file and the user
+provided configuration file. The final configuration is obtained by
+overriding the options by the values provided by the user.
+
+Users can provide additional configuration files, by specifying the
+location of their file in the environmental variable
+`SITS_CONFIG_USER_FILE` or as parameter to this function.
+To see the key entries and contents of the current configuration
+values, use `sits_config_show()`.
+
+Args:
+    config_user_file (str | pathlib.Path): YAML user configuration file
+        (a file with a "yml" extension).
+
+Returns:
+    SITStructureData: Called for side effects.
+
+Examples:
+    from pysits import *
+
+    yaml_user_file = r_package_dir("extdata/config_user_example.yml", package="sits")
+    sits_config(config_user_file=yaml_user_file)
diff --git a/pysits/docs/content/sits_config_show.md b/pysits/docs/content/sits_config_show.md
new file mode 100644
index 0000000..bd0d1dc
--- /dev/null
+++ b/pysits/docs/content/sits_config_show.md
@@ -0,0 +1,13 @@
+Show current sits configuration
+
+Prints the current sits configuration options. To show specific
+configuration options for a source, a collection, or a palette, users
+can inform the corresponding keys to `source` and `collection`.
+
+Returns:
+    None: called for side effects.
+
+Examples:
+    from pysits import *
+
+    sits_config_show()
diff --git a/pysits/docs/content/sits_config_user_file.md b/pysits/docs/content/sits_config_user_file.md
new file mode 100644
index 0000000..6708937
--- /dev/null
+++ b/pysits/docs/content/sits_config_user_file.md
@@ -0,0 +1,19 @@
+Create a user configuration file.
+
+Creates a user configuration file.
+
+Args:
+    file_path (str | pathlib.Path): file to store the user
+        configuration file.
+    overwrite (bool): replace current configuration file?
+
+Returns:
+    None: called for side effects.
+
+Examples:
+    from pysits import *
+    import tempfile
+    import os
+
+    user_file = os.path.join(tempfile.gettempdir(), "my_config_file.yml")
+    sits_config_user_file(user_file)
diff --git a/pysits/docs/content/sits_contrastive_learning.md b/pysits/docs/content/sits_contrastive_learning.md
new file mode 100644
index 0000000..6eb1b5b
--- /dev/null
+++ b/pysits/docs/content/sits_contrastive_learning.md
@@ -0,0 +1,84 @@
+Supervised contrastive learning for time series
+
+Supervised contrastive (SupCon) pre-training of encoders, based on the
+loss proposed by Khosla et al. (2020). The method learns an embedding space in
+which time-series samples sharing the same label are pulled together while
+samples with different labels are pushed apart.
+For each batch, two views per sample are created by pairing every anchor with a
+same-class sample. Both views are passed through a shared encoder and
+projection head and are L2-normalized. The resulting `2B` embeddings are pooled
+and each embedding is contrasted against every other embedding (excluding
+itself) using temperature-scaled cosine similarity, with positives defined by
+matching labels.
+After pre-training, the projection head is discarded and only the encoder is
+kept for downstream use via `sits_encode`.
+The function can be used in two ways:
+- If `samples` is provided, it trains immediately and returns an encoder-ready
+  model object (see Value).
+- If `samples = None`, it returns a training function with signature
+  `function(samples)` that can be passed to `sits_pre_train` or called later.
+
+The supervised contrastive (SupCon) loss generalizes the InfoNCE/NT-Xent
+objective by allowing multiple positives per anchor. For each of the `2B`
+L2-normalized embeddings in a batch, all other embeddings that share the
+anchor's label act as positives and the remaining embeddings act as negatives.
+Pairwise cosine similarities are divided by the temperature `scaling` and
+combined in a log-softmax form; the loss is averaged over the positives of each
+anchor and then over the batch.
+The `scaling` parameter (temperature) controls the sharpness of the similarity
+distribution. Lower temperatures sharpen it and place more weight on the
+hardest negatives, encouraging stronger separation at the cost of noisier
+gradients; higher temperatures soften it. The default of `0.07` follows Khosla
+et al. (2020). Larger batches also help, since they expose more positives and
+negatives per anchor and yield a stronger contrastive signal.
+
+Args:
+    samples (SITSTimeSeriesModel): Samples object. If `None` (default),
+        returns a training function. If provided, triggers immediate
+        training. Base data samples are not supported.
+    embedding_dim (int): Dimensionality of the encoder embedding (exported
+        features). Default: 64.
+    proj_dim (int): Dimensionality of the projection head output used only
+        during pre-training (discarded afterwards). Default: 128.
+    scaling (float): Scaling for the contrastive loss. Lower values sharpen
+        the similarity distribution. Default: 0.07.
+    num_pairs (int | None): Total number of pairs to form. When `None`
+        (default), one pair is formed per sample.
+    encoder_model (SITSMachineLearningMethod): Deep learning method that
+        takes time series as input and produces latent representations that
+        are used to compute the loss function (suggested options:
+        `sits_tempcnn()`, `sits_lighttae()`, `sits_resnet()`). Default:
+        `sits_tempcnn()`.
+    epochs (int): Maximum number of training epochs.
+    batch_size (int): Batch size for training. Larger batches provide more
+        positives/negatives per sample. Default: 128.
+    validation_split (float): Fraction of samples held out for validation
+        loss monitoring, in the range (0, 1).
+    optimizer: A `torch` optimizer constructor (default:
+        `torch::optim_adamw`).
+    opt_hparams (dict): Optimizer hyperparameters. Common entries: `lr`,
+        `eps`, `weight_decay`.
+    lr_decay_epochs (int): Step size (in epochs) for LR decay.
+    lr_decay_rate (float): Multiplicative LR decay factor.
+    patience (int): Early-stopping patience (epochs without improvement).
+    min_delta (float): Minimum improvement required to reset the patience
+        counter.
+    verbose (bool): Print training progress?
+    seed (int): Random seed for reproducibility.
+
+Returns:
+    R: If `samples = None`, a training function with signature
+    `function(samples)` that trains a supervised contrastive model and returns
+    a pretrained encoder. If `samples` is provided, the result of applying the
+    training function to `samples` directly.
+
+Examples:
+    from pysits import *
+
+    model = sits_pre_train(
+        samples_modis_ndvi,
+        sits_contrastive_learning(
+            embedding_dim=32,
+            epochs=20
+        )
+    )
diff --git a/pysits/docs/content/sits_cube.md b/pysits/docs/content/sits_cube.md
new file mode 100644
index 0000000..dd88222
--- /dev/null
+++ b/pysits/docs/content/sits_cube.md
@@ -0,0 +1,97 @@
+Create sits data cubes from local files or cloud-based image collections.
+
+Builds a data cube -- a `pandas.DataFrame` describing spatial and temporal
+image data -- from one of several sources. This single function handles four
+distinct cases, selected by the combination of arguments you provide:
+
+- Local raster cube: images already downloaded from a known cloud collection
+  or created by `sits`, read from a local directory (`data_dir`).
+- STAC cube: image collections accessible through the STAC protocol, selected
+  by spatial (`roi`, `tiles`) and temporal (`start_date`, `end_date`)
+  restrictions.
+- Results cube: local files produced by `sits` operations that generate
+  results (for example probability cubes and class cubes).
+- Vector cube: local files that include a vector file produced by a
+  segmentation algorithm, merged with a raster cube.
+
+Args:
+    source (str): Data source: one of `"AWS"`, `"BDC"`, `"CDSE"`,
+            `"DEAFRICA"`, `"DEAUSTRALIA"`, `"HLS"`, `"PLANETSCOPE"`,
+            `"MPC"`, `"SDC"` or `"USGS"`. For local, results and vector
+            cubes this is the source from which the original data was
+            downloaded.
+    collection (str): Image collection in the data source. Use
+            `sits_list_collections()` to find supported collections.
+    bands (list[str]): Spectral bands and indices to include in the cube
+            (optional). For a results cube these are results bands to be
+            retrieved (`"probs"`, `"bayes"`, `"variance"`, `"class"`,
+            `"uncertainty"`). Use `sits_list_collections()` to find the
+            bands available for each STAC collection.
+    tiles (list[str]): Tiles from the collection to include in the cube.
+    roi (dict): Region of interest for STAC cubes.
+    crs (str): The Coordinate Reference System (CRS) of the `roi`.
+    start_date (str): Initial date to include images from the collection in
+            the cube, in `YYYY-MM-DD` format (optional).
+    end_date (str): Final date to include images from the collection in the
+            cube, in `YYYY-MM-DD` format (optional).
+    orbit (str): Orbit name (`"ascending"`, `"descending"`) for SAR cubes.
+    platform (str): Optional parameter specifying the platform for the
+            `"LANDSAT"` collection. Options: `Landsat-5`, `Landsat-7`,
+            `Landsat-8`, `Landsat-9`.
+    data_dir (str | pathlib.Path): Local directory where images are stored
+            (local, results and vector cubes).
+    labels (dict): Labels associated to the classes (results cube).
+    raster_cube (SITSCubeModel): Raster cube to be merged with vector data
+            (vector cube).
+    vector_dir (str | pathlib.Path): Local directory where vector files are
+            stored (vector cube).
+    vector_band (str): Band for the vector cube (`"segments"`, `"probs"`,
+            `"class"`). Deprecated and will be removed in future versions;
+            the vector data cube type is now defined from the `raster_cube`
+            object.
+    parse_info (list[str]): Parsing information for local files.
+    version (str): Version of the classified and/or labelled files (results
+            and vector cubes).
+    delim (str): Delimiter for parsing local files (default `"_"`).
+    multicores (int): Number of workers for parallel processing
+            (min = 1, max = 2048).
+    memsize (int): Memory available in GB (results cube).
+    progress (bool): Whether to show a progress bar.
+    **kwargs (dict): Other parameters passed to specific cube types.
+
+Returns:
+    SITSCubeModel: A data cube describing the contents of the images.
+
+Notes:
+    The specific cube type is inferred from the combination of arguments.
+    Local raster cubes use `data_dir`; STAC cubes use spatial and temporal
+    restrictions such as `roi`, `tiles`, `start_date` and `end_date`;
+    results cubes use `data_dir` together with results `bands` and
+    `labels`; vector cubes use `raster_cube` together with `vector_dir`.
+
+Examples:
+    from pysits import *
+
+    # Create a local data cube from MODIS files bundled with sits
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # Inspect the cube
+    print(sits_bands(cube))
+    print(sits_timeline(cube))
+    plot(cube)
+
+    # Create a STAC-based data cube restricted by tiles and dates
+    s2_cube = sits_cube(
+        source="MPC",
+        collection="SENTINEL-2-L2A",
+        tiles="20LKP",
+        bands=["B05", "CLOUD"],
+        start_date="2018-07-18",
+        end_date="2018-08-23"
+    )
+    print(sits_bands(s2_cube))
diff --git a/pysits/docs/content/sits_cube_copy.md b/pysits/docs/content/sits_cube_copy.md
new file mode 100644
index 0000000..8e94425
--- /dev/null
+++ b/pysits/docs/content/sits_cube_copy.md
@@ -0,0 +1,85 @@
+Copy the images of a cube to a local directory
+
+This function downloads the images of a cube in parallel. A region of
+interest (`roi`) can be provided to crop the images and a resolution
+(`res`) to resample the bands. `sits_cube_copy` is useful to improve
+processing time in the regularization operation.
+
+Args:
+    cube (SITSCubeModel): A data cube.
+    roi (str | pathlib.Path | geopandas.GeoDataFrame | dict): Region of
+        interest. Either: 1. A path to a shapefile with polygons; 2. A
+        `geopandas.GeoDataFrame`; 3. A `dict` (`"lon_min"`, `"lat_min"`,
+        `"lon_max"`, `"lat_max"`) in WGS84; 4. A `dict` (`"xmin"`,
+        `"xmax"`, `"ymin"`, `"ymax"`) with XY coordinates in the
+        projection of the input cube.
+    res (int): Output spatial resolution of the images. Default is None.
+    crs (str): The Coordinate Reference System (CRS) of the roi. (see
+        details below).
+    n_tries (int): Number of attempts to download the same image. Default
+        is 3.
+    multicores (int): Number of cores for parallel downloading (min = 1,
+        max = 2048).
+    output_dir (str | pathlib.Path): Output directory where images will be
+        saved.
+    progress (bool): Show progress bar?
+    **kwargs (dict): Additional arguments.
+
+Returns:
+    SITSCubeModel: Copy of input data cube.
+    The main `sits` classification workflow has the following steps:
+    1. `sits_cube`: selects a ARD image collection from a cloud provider.
+    2. `sits_cube_copy`: copies an ARD image collection from a cloud
+       provider to a local directory for faster processing.
+    3. `sits_regularize`: create a regular data cube from an ARD image
+       collection.
+    4. `sits_apply`: create new indices by combining bands of a regular
+       data cube (optional).
+    5. `sits_get_data`: extract time series from a regular data cube based
+       on user-provided labelled samples.
+    6. `sits_train`: train a machine learning model based on image time
+       series.
+    7. `sits_classify`: classify a data cube using a machine learning
+       model and obtain a probability cube.
+    8. `sits_smooth`: post-process a probability cube using a spatial
+       smoother to remove outliers and increase spatial consistency.
+    9. `sits_label_classification`: produce a classified map by selecting
+       the label with the highest probability from a smoothed cube.
+    The `roi` parameter is used to crop cube images. To define a `roi`
+    use one of:
+    - A path to a shapefile with polygons;
+    - A `geopandas.GeoDataFrame`;
+    - A spatial extent object;
+    - A `dict` (`"lon_min"`, `"lat_min"`, `"lon_max"`, `"lat_max"`) in
+      WGS84;
+    - A `dict` (`"xmin"`, `"xmax"`, `"ymin"`, `"ymax"`) with XY
+      coordinates.
+    Defining a region of interest using a spatial extent or XY values not
+    in WGS84 requires the `crs` parameter to be specified.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # Creating a sits cube from BDC
+    bdc_cube = sits_cube(
+        source="BDC",
+        collection="CBERS-WFI-16D",
+        tiles=["007004", "007005"],
+        bands=["B15", "CLOUD"],
+        start_date="2018-01-01",
+        end_date="2018-01-12"
+    )
+    # Downloading images to a temporary directory
+    cube_local = sits_cube_copy(
+        cube=bdc_cube,
+        output_dir=tempfile.gettempdir(),
+        roi={
+            "lon_min": -46.5,
+            "lat_min": -45.5,
+            "lon_max": -15.5,
+            "lat_max": -14.6
+        },
+        multicores=2,
+        res=250
+    )
diff --git a/pysits/docs/content/sits_encode.md b/pysits/docs/content/sits_encode.md
new file mode 100644
index 0000000..6ab1c6c
--- /dev/null
+++ b/pysits/docs/content/sits_encode.md
@@ -0,0 +1,57 @@
+Encode data using a pre-trained deep learning encoder.
+
+Encodes either a regular raster data cube or a set of time series using a
+pre-trained deep learning encoder returned by `sits_pre_train`. The
+behavior depends on the type of `data` provided.
+
+When `data` is a raster data cube, the output is an embeddings cube with
+the same tiling as the input cube. Each tile is written as a multiband
+raster, where each band corresponds to one embedding dimension produced
+by the encoder.
+
+When `data` is a set of time series, the output preserves the input
+structure and replaces the `time_series` column with the corresponding
+embeddings.
+
+Args:
+    data (SITSCubeModel | SITSTimeSeriesModel): Data to be encoded. Either
+        a regular raster data cube or a set of time series.
+    encoder (SITSRepresentationLearningMethod): Encoder returned by
+        `sits_pre_train`.
+    roi (str | pathlib.Path | geopandas.GeoDataFrame | dict): Optional
+        region of interest used to restrict processing (raster cube only).
+        It may be provided as: (1) a path to a polygon shapefile; (2) a
+        `geopandas.GeoDataFrame` with `POLYGON` or `MULTIPOLYGON`
+        geometry; (3) a named bounding box in WGS84 with `xmin`, `xmax`,
+        `ymin`, `ymax`; or (4) a named lon/lat bounding box with
+        `lon_min`, `lon_max`, `lat_min`, `lat_max`.
+    impute_fn: Imputation function used to interpolate missing values in
+        each pixel time series (default: `impute_linear`).
+    start_date (str): Optional start date for temporal filtering
+        (YYYY-MM-DD), raster cube only. Defaults to the cube start date.
+    end_date (str): Optional end date for temporal filtering (YYYY-MM-DD),
+        raster cube only. Defaults to the cube end date.
+    memsize (int): Memory available for processing in GB (minimum 1),
+        raster cube only.
+    multicores (int): Number of CPU cores used for processing (minimum 1;
+        maximum 2048 for time series).
+    gpu_memory (int): GPU memory available for encoding in GB (minimum 1).
+    batch_size (int): Batch size used when encoding on GPU.
+    block_size (dict): Size of the block read and written by each worker
+        (raster cube only). A named vector with `nrows` and `ncols`.
+        Default is `None`, which computes an optimal block size from
+        `memsize`, `multicores` and the internal block size of the raster
+        files.
+    output_dir (str | pathlib.Path): Directory where output files will be
+        written (raster cube only).
+    verbose (bool): If `True`, print processing time information (raster
+        cube only).
+    progress (bool): If `True`, show a progress bar.
+    **kwargs (dict): Additional arguments passed to lower-level encoding
+        routines.
+
+Returns:
+    SITSTimeSeriesModel: For a raster cube, an embeddings cube written to
+    `output_dir`, with one multiband raster per tile date. For a set of
+    time series, a `SITSTimeSeriesModel` with `time_series` containing the
+    embeddings produced by `encoder`.
diff --git a/pysits/docs/content/sits_formula_linear.md b/pysits/docs/content/sits_formula_linear.md
new file mode 100644
index 0000000..281acd5
--- /dev/null
+++ b/pysits/docs/content/sits_formula_linear.md
@@ -0,0 +1,30 @@
+Define a linear formula for classification models
+
+Provides a symbolic description of a fitting model. Tells the model to do a
+linear transformation of the input values. The `predictors_index` parameter
+informs the positions of fields corresponding to formula independent variables.
+If no value is given, that all fields will be used as predictors.
+
+Args:
+    predictors_index (list[int]): Index of the valid columns whose names are
+        used to compose formula (default: -2:0).
+
+Returns:
+    R: A function that computes a valid formula using a linear function.
+
+Examples:
+    from pysits import *
+
+    # Example of training a model for time series classification
+    # Retrieve the samples for Mato Grosso
+    # train an SVM model
+    ml_model = sits_train(samples_modis_ndvi,
+        ml_method=sits_svm(formula=sits_formula_logref())
+    )
+    # classify the point
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    # classify the point
+    point_class = sits_classify(
+        data=point_ndvi, ml_model=ml_model
+    )
+    plot(point_class)
diff --git a/pysits/docs/content/sits_formula_logref.md b/pysits/docs/content/sits_formula_logref.md
new file mode 100644
index 0000000..e2dba76
--- /dev/null
+++ b/pysits/docs/content/sits_formula_logref.md
@@ -0,0 +1,32 @@
+Define a loglinear formula for classification models
+
+A function to be used as a symbolic description of some fitting models such
+as svm and random forest. This function tells the models to do a log
+transformation of the inputs. The `predictors_index` parameter informs the
+positions of `tb` fields corresponding to formula independent variables. If
+no value is given, the default is `None`, a value indicating that all fields
+will be used as predictors.
+
+Args:
+    predictors_index (list[int]): Index of the valid columns to compose
+        formula (default: -2:0).
+
+Returns:
+    R: A function that computes a valid formula using a log function.
+
+Examples:
+    from pysits import *
+
+    # Example of training a model for time series classification
+    # Retrieve the samples for Mato Grosso
+    # train an SVM model
+    ml_model = sits_train(samples_modis_ndvi,
+        ml_method=sits_svm(formula=sits_formula_logref())
+    )
+    # classify the point
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    # classify the point
+    point_class = sits_classify(
+        data=point_ndvi, ml_model=ml_model
+    )
+    plot(point_class)
diff --git a/pysits/docs/content/sits_geo_dist.md b/pysits/docs/content/sits_geo_dist.md
new file mode 100644
index 0000000..7d5a88c
--- /dev/null
+++ b/pysits/docs/content/sits_geo_dist.md
@@ -0,0 +1,42 @@
+Compute the minimum distances among samples and prediction points.
+
+Compute the minimum distances among samples and samples to prediction points,
+following the approach proposed by Meyer and Pebesma(2022).
+
+Args:
+    samples (SITSTimeSeriesModel): Time series.
+    roi (str | pathlib.Path | geopandas.GeoDataFrame): A region of interest
+        (ROI), either a file containing a shapefile or a
+        `geopandas.GeoDataFrame`.
+    n (int): Maximum number of samples to consider.
+    crs (str): CRS of the `samples`.
+
+Returns:
+    SITSFrame: sample-to-sample and sample-to-prediction distances.
+
+Notes:
+    As pointed out by Meyer and Pebesma, many classifications using machine
+    learning assume that the reference data are independent and well-
+    distributed in space. In practice, many training samples are strongly
+    concentrated in some areas, and many large areas have no samples. This
+    function compares two distributions:
+    1. The distribution of the spatial distances of reference data to their
+       nearest neighbor (sample-to-sample.
+    2. The distribution of distances from all points of study area to the
+       nearest reference data point (sample-to-prediction).
+
+Examples:
+    from pysits import *
+    import geopandas as gpd
+
+    # read a shapefile for the state of Mato Grosso, Brazil
+    mt_shp = r_package_dir("extdata/shapefiles/mato_grosso/mt.shp", package="sits")
+    # convert to a geopandas object
+    mt_sf = gpd.read_file(mt_shp)
+    # calculate sample-to-sample and sample-to-prediction distances
+    distances = sits_geo_dist(
+        samples=samples_modis_ndvi,
+        roi=mt_sf
+    )
+    # plot sample-to-sample and sample-to-prediction distances
+    plot(distances)
diff --git a/pysits/docs/content/sits_get_class.md b/pysits/docs/content/sits_get_class.md
new file mode 100644
index 0000000..fd50d13
--- /dev/null
+++ b/pysits/docs/content/sits_get_class.md
@@ -0,0 +1,56 @@
+Get values from classified maps
+
+Given a set of lat/long locations and a classified cube, retrieve the class of
+each point. This function is useful to obtain values from classified cubes for
+accuracy estimates.
+
+Args:
+    cube (SITSCubeModel): Classified data cube.
+    samples (SITSTimeSeriesModel | geopandas.GeoDataFrame | str | pathlib.Path | pandas.DataFrame): Location of the
+        samples to be retrieved. Either a `SITSTimeSeriesModel`, a
+        `geopandas.GeoDataFrame`, the name of a shapefile or CSV file, or a
+        `pandas.DataFrame` with columns "longitude" and "latitude".
+
+Returns:
+    SITSFrame: With columns <longitude, latitude, start_date, end_date,
+        label>.
+
+Notes:
+    There are four ways of specifying data to be retrieved using the `samples`
+    parameter: (a) CSV file: a CSV file with columns `longitude`, `latitude`;
+    (b) SHP file: a shapefile in POINT geometry; (c) a `SITSTimeSeriesModel`;
+    (d) a `geopandas.GeoDataFrame` with POINT or geometry; (e) a
+    `pandas.DataFrame` with `longitude` and `latitude`.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # obtain the a set of points for sampling
+    ground_truth = r_package_dir("extdata/samples/samples_sinop_crop.csv", package="sits")
+    # get the classification values for a selected set of locations
+    labels_samples = sits_get_class(label_cube, ground_truth)
diff --git a/pysits/docs/content/sits_get_data.md b/pysits/docs/content/sits_get_data.md
index c0cfcbb..7ce84b8 100644
--- a/pysits/docs/content/sits_get_data.md
+++ b/pysits/docs/content/sits_get_data.md
@@ -1,85 +1,80 @@
-Get time series from data cubes and cloud services
-
-Description:
-    Retrieve a set of time series from a data cube and and put the result 
-    in a sits tibble, which contains both the satellite image time series 
-    and their metadata.
-
-    There are five options for specifying the input samples parameter:
-
-    A CSV file: see sits_get_data.csv.
-
-    A shapefile: see sits_get_data.shp.
-
-    An sf object: see sits_get_data.sf.
-
-    A sits tibble: see sits_get_data.sits.
-
-    A data.frame: see see sits_get_data.data.frame.
-
-Usage:
-    sits_get_data(cube, samples, ...)
-
-    ## Default S3 method:
-    sits_get_data(cube, samples, ...)
+Get time series from a data cube.
+
+Retrieve a set of time series from a data cube and put the result in a
+`SITSTimeSeriesModel`, which contains both the satellite image time
+series and their metadata. The samples to be retrieved can be provided
+in several forms, and the accepted parameters vary slightly depending
+on the type of `samples`.
+
+The `samples` parameter may be one of the following:
+
+- A path to a CSV file (extension ".csv") with mandatory columns
+  `longitude`, `latitude`, `label`, `start_date` and `end_date`.
+- A `pandas.DataFrame` with mandatory columns `longitude` and
+  `latitude`, and optional columns `start_date`, `end_date` and
+  `label`.
+- A `geopandas.GeoDataFrame` in POINT or POLYGON geometry.
+- A path to a shapefile (extension ".shp") that is a valid POINT or
+  POLYGON shapefile.
+- A valid `SITSTimeSeriesModel` with columns `longitude`, `latitude`,
+  `start_date`, `end_date` and `label`.
+
+For spatial inputs (`geopandas.GeoDataFrame` objects and shapefiles),
+if `start_date` and `end_date` are not informed, the function uses
+these dates from the cube.
 
 Args:
-    cube: 
-        Data cube from where data is to be retrieved. 
-        (tibble of class "raster_cube").
-
-    samples: 
-        Location of the samples to be retrieved. Either a tibble of class 
-        "sits", an "sf" object, the name of a shapefile or csv file, or a 
-        data.frame with columns "longitude" and "latitude".
-
-    **kwargs: 
-        Specific parameters for each input.
+    cube (SITSCubeModel): Data cube from where data is to be retrieved.
+    samples (SITSTimeSeriesModel | geopandas.GeoDataFrame | pandas.DataFrame | str | pathlib.Path): Location of the samples
+        to be retrieved. Either a `SITSTimeSeriesModel`, a
+        `geopandas.GeoDataFrame`, the name of a shapefile or csv file,
+        or a `pandas.DataFrame` with columns "longitude" and
+        "latitude".
+    bands (list[str]): Bands to be retrieved - optional.
+    start_date (str): Start of the interval for the time series -
+        optional (date in "YYYY-MM-DD" format). Applies to
+        `pandas.DataFrame`, `geopandas.GeoDataFrame` and shapefile
+        inputs.
+    end_date (str): End of the interval for the time series - optional
+        (date in "YYYY-MM-DD" format). Applies to `pandas.DataFrame`,
+        `geopandas.GeoDataFrame` and shapefile inputs.
+    label (str): Label to be assigned to all time series if a `label`
+        column is not provided in the input.
+    label_attr (str): Attribute in the `geopandas.GeoDataFrame` or
+        shapefile to be used as a polygon label.
+    n_sam_pol (int): Number of samples per polygon to be read for
+        POLYGON or MULTIPOLYGON inputs.
+    pol_avg (bool): Summarize samples for each polygon?
+    sampling_type (str): Spatial sampling type: random, hexagonal,
+        regular, or Fibonacci.
+    crs (str): The samples CRS. Default is "EPSG:4326".
+    impute_fn: Imputation function to remove NA.
+    multicores (int): Number of threads to process the time series
+        (with min = 1 and max = 2048).
+    progress (bool): Show progress bar?
+    **kwargs (dict): Specific parameters for each kind of input.
 
 Returns:
-    A tibble of class "sits" with set of time series 
-    <longitude, latitude, start_date, end_date, label, time_series>.
-
-Note:
-    The main sits classification workflow has the following steps:
-
-    sits_cube: selects a ARD image collection from a cloud provider.
-
-    sits_cube_copy: copies an ARD image collection from a cloud provider 
-    to a local directory for faster processing.
-
-    sits_regularize: create a regular data cube from an ARD image collection.
-
-    sits_apply: create new indices by combining bands of a regular data 
-    cube (optional).
-
-    sits_get_data: extract time series from a regular data cube based on 
-    user-provided labelled samples.
-
-    sits_train: train a machine learning model based on image time series.
-
-    sits_classify: classify a data cube using a machine learning model and 
-    obtain a probability cube.
-
-    sits_smooth: post-process a probability cube using a spatial smoother 
-    to remove outliers and increase spatial consistency.
-
-    sits_label_classification: produce a classified map by selecting the 
-    label with the highest probability from a smoothed cube.
-
-    To be able to build a machine learning model to classify a data cube, 
-    one needs to use a set of labelled time series. These time series are 
-    created by taking a set of known samples, expressed as labelled points 
-    or polygons. This sits_get_data function uses these samples to extract 
-    time series from a data cube. It needs a cube parameter which points to 
-    a regularized data cube, and a samples parameter that describes the 
-    locations of the training set.
-
-Author(s):
-    Felipe Carlos, efelipecarlos@gmail.com
-
-    Felipe Carvalho, felipe.carvalho@inpe.br
-
-    Gilberto Camara, gilberto.camara@inpe.br
-
-    Rolf Simoes, rolfsimoes@gmail.com
+    SITSTimeSeriesModel: The set of time series and metadata:
+        <longitude, latitude, start_date, end_date, label,
+        time_series>.
+
+Examples:
+    from pysits import *
+    import pandas as pd
+
+    # reading a lat/long from a local cube
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    raster_cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # obtain a set of samples defined by a lat/long point
+    samples = pd.DataFrame({
+        "longitude": [-55.66738],
+        "latitude": [-11.76990]
+    })
+    points = sits_get_data(cube=raster_cube, samples=samples)
+    print(points)
diff --git a/pysits/docs/content/sits_get_probs.md b/pysits/docs/content/sits_get_probs.md
new file mode 100644
index 0000000..c7ab5a9
--- /dev/null
+++ b/pysits/docs/content/sits_get_probs.md
@@ -0,0 +1,53 @@
+Get values from probability maps
+
+Given a set of lat/long locations and a probability cube, retrieve the prob
+values of each point. This function is useful to estimate probability
+distributions and to assess the differences between classifiers.
+
+Args:
+    cube (SITSCubeModel): Probability data cube.
+    samples (SITSTimeSeriesModel | geopandas.GeoDataFrame | str | pathlib.Path | pandas.DataFrame):
+        Location of the samples to be retrieved. Either a
+        `SITSTimeSeriesModel`, a `geopandas.GeoDataFrame` with POINT
+        geometry, the location of a POINT shapefile, the location of a
+        CSV file with columns `longitude` and `latitude`, or a
+        `pandas.DataFrame` with columns `longitude` and `latitude`.
+    window_size (int): Size of window around pixel (optional).
+    **kwargs (dict): Additional arguments.
+
+Returns:
+    SITSFrameNested: A table with columns <longitude, latitude, values> in
+        case no windows are requested and <longitude, latitude, neighbors>
+        in case windows are requested.
+
+Notes:
+    There are four ways of specifying data to be retrieved using the `samples`
+    parameter:
+    - CSV: a CSV file with columns `longitude`, `latitude`.
+    - SHP: a shapefile in POINT geometry.
+    - `geopandas.GeoDataFrame`: an object with POINT geometry.
+    - `SITSTimeSeriesModel`: a valid set of time series.
+    - `pandas.DataFrame`: a data frame with `longitude` and `latitude`.
+
+Examples:
+    from pysits import *
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempdir()
+    )
+    # obtain the a set of points for sampling
+    ground_truth = r_package_dir("extdata/samples/samples_sinop_crop.csv",
+        package="sits"
+    )
+    # get the classification values for a selected set of locations
+    probs_samples = sits_get_probs(probs_cube, ground_truth)
diff --git a/pysits/docs/content/sits_kfold_validate.md b/pysits/docs/content/sits_kfold_validate.md
new file mode 100644
index 0000000..f0e515b
--- /dev/null
+++ b/pysits/docs/content/sits_kfold_validate.md
@@ -0,0 +1,58 @@
+Cross-validate time series samples
+
+Splits the set of time series into training and validation and perform k-fold
+cross-validation.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series.
+    folds (int): Number of partitions to create.
+    ml_method (SITSMachineLearningMethod): Machine learning method.
+    impute_fn: Imputation function to remove NA.
+    multicores (int): Number of cores to process in parallel.
+    gpu_memory (int): Memory available in GPU in GB (default = 4).
+    batch_size (int): Batch size for GPU classification.
+    progress (bool): Show progress bar?
+    **kwargs (dict): Additional arguments.
+
+Returns:
+    resolve_and_invoke_accuracy_class: A confusion matrix object to be used
+        for validation assessment.
+
+Notes:
+    Cross-validation is a technique for assessing how the results of a
+    statistical analysis will generalize to an independent data set. It is
+    mainly used in settings where the goal is prediction, and one wants to
+    estimate how accurately a predictive model will perform. One round of
+    cross-validation involves partitioning a sample of data into complementary
+    subsets, performing the analysis on one subset (called the training set),
+    and validating the analysis on the other subset (called the validation set
+    or testing set).
+    The k-fold cross validation method involves splitting the dataset into
+    k-subsets. For each subset is held out while the model is trained on all
+    other subsets. This process is completed until accuracy is determine for
+    each instance in the dataset, and an overall accuracy estimate is provided.
+    This function returns the confusion matrix, and Kappa values.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # A dataset containing a tibble with time series samples
+    # for the Mato Grosso state in Brasil
+    # create a list to store the results
+    results = []
+    # accuracy assessment random forest
+    acc_rfor = sits_kfold_validate(
+        samples_modis_ndvi,
+        folds=5,
+        ml_method=sits_rfor()
+    )
+    # use a name
+    acc_rfor.name = "Rfor"
+    # put the result in a list
+    results.append(acc_rfor)
+    # save to xlsx file
+    sits_to_xlsx(
+        results,
+        file=tempfile.NamedTemporaryFile(prefix="accuracy_mato_grosso_dl_", suffix=".xlsx").name
+    )
diff --git a/pysits/docs/content/sits_label_classification.md b/pysits/docs/content/sits_label_classification.md
new file mode 100644
index 0000000..af94f53
--- /dev/null
+++ b/pysits/docs/content/sits_label_classification.md
@@ -0,0 +1,91 @@
+Build a labelled image from a probability cube
+
+Takes a set of classified raster layers with probabilities, and labels them
+based on the maximum probability for each pixel. This function is the final
+step of the main land classification workflow.
+When the input is a segmented probability cube (produced by `sits_classify`
+from a segmented cube), this function performs segment-based labeling:
+pixel-level probabilities are aggregated inside each segment using the method
+specified by `label_method`. All pixels within a segment receive the same
+class label in the output raster. A GPKG file with segment summaries
+(including a `class` column) is written automatically.
+
+Args:
+    cube (SITSCubeModel): Classified probability data cube (pixel-based or
+        segment-based).
+    memsize (int): Maximum overall memory (in GB) to label the
+        classification.
+    multicores (int): Number of workers to label the classification in
+        parallel.
+    output_dir (str | pathlib.Path): Output directory for classified files.
+    version (str): Version of resulting image (in the case of multiple
+        runs).
+    progress (bool): Show progress bar?
+    label_method (str): Decision method for segment-based labeling. One of
+        "mean" (default), "median", or "majority". Only used when input is a
+        segment-based probability cube.
+    **kwargs (dict): Configuration parameters for exact extraction of
+        segment values.
+
+Returns:
+    SITSCubeModel: A data cube with an image with the classified map. When
+        input is a segment-based probability cube, the output preserves
+        vector information.
+
+Notes:
+    The main `sits` classification workflow has the following steps:
+    1. `sits_cube`: selects a ARD image collection from a cloud provider.
+    2. `sits_cube_copy`: copies an ARD image collection from a cloud provider
+       to a local directory for faster processing.
+    3. `sits_regularize`: create a regular data cube from an ARD image
+       collection.
+    4. `sits_apply`: create new indices by combining bands of a regular data
+       cube (optional).
+    5. `sits_get_data`: extract time series from a regular data cube based on
+       user-provided labelled samples.
+    6. `sits_train`: train a machine learning model based on image time series.
+    7. `sits_classify`: classify a data cube using a machine learning model and
+       obtain a probability cube.
+    8. `sits_smooth`: post-process a probability cube using a spatial smoother
+       to remove outliers and increase spatial consistency.
+    9. `sits_label_classification`: produce a classified map by selecting the
+       label with the highest probability from a smoothed cube.
+    The OBIA workflow adds segmentation before classification:
+    1. `sits_segment`: segment the raster cube to produce a vector_cube.
+    2. `sits_classify`: classify pixel-level probabilities, preserving vector
+       support.
+    3. `sits_label_classification`: aggregate probabilities per segment and
+       assign class labels.
+    Please refer to the sits documentation available in
+    https://e-sensing.github.io/sitsbook/ for detailed examples.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_labels.md b/pysits/docs/content/sits_labels.md
new file mode 100644
index 0000000..0c3f6ce
--- /dev/null
+++ b/pysits/docs/content/sits_labels.md
@@ -0,0 +1,35 @@
+Get labels associated to a data set
+
+Finds labels in a time series set or data cube
+
+Args:
+    data (SITSTimeSeriesModel | SITSTimeSeriesPatternsModel | SITSCubeModel | SITSMachineLearningMethod): time series, patterns, data cube, or trained model.
+
+Returns:
+    list: the labels of the input data.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # get the labels for a time series set
+    labels_ts = sits_labels(samples_modis_ndvi)
+    # get labels for a set of patterns
+    labels_pat = sits_labels(sits_patterns(samples_modis_ndvi))
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # get labels for the model
+    labels_mod = sits_labels(rfor_model)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # get the labels for a probs cube
+    labels_probs = sits_labels(probs_cube)
diff --git a/pysits/docs/content/sits_labels_summary.md b/pysits/docs/content/sits_labels_summary.md
new file mode 100644
index 0000000..3ceec10
--- /dev/null
+++ b/pysits/docs/content/sits_labels_summary.md
@@ -0,0 +1,16 @@
+Inform label distribution of a set of time series
+
+Describes labels in a set of time series.
+
+Args:
+    data (SITSTimeSeriesModel): Valid time series.
+
+Returns:
+    SITSFrame: A table with the frequency of each label.
+
+Examples:
+    from pysits import *
+
+    # read a tibble with 400 samples of Cerrado and 346 samples of Pasture
+    # print the labels
+    sits_labels_summary(cerrado_2classes)
diff --git a/pysits/docs/content/sits_lightgbm.md b/pysits/docs/content/sits_lightgbm.md
new file mode 100644
index 0000000..2370aef
--- /dev/null
+++ b/pysits/docs/content/sits_lightgbm.md
@@ -0,0 +1,42 @@
+Train light gradient boosting model
+
+Use LightGBM algorithm to classify samples. This function is a front-end to
+the `lightgbm` package. LightGBM (short for Light Gradient Boosting Machine)
+is a gradient boosting framework developed by Microsoft that's designed for
+fast, scalable, and efficient training of decision tree-based models. It is
+widely used in machine learning for classification, regression, ranking, and
+other tasks, especially with large-scale data.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    boosting_type (str): Type of boosting algorithm (default = "gbdt").
+    objective (str): Aim of the classifier (default = "multiclass").
+    min_samples_leaf (int): Minimal number of data in one leaf. Can be used
+        to deal with over-fitting.
+    max_depth (int): Limit the max depth for tree model.
+    learning_rate (float): Shrinkage rate for leaf-based algorithm.
+    num_iterations (int): Number of iterations to train the model.
+    n_iter_no_change (int): Number of iterations without improvements until
+        training stops.
+    validation_split (float): Fraction of the training data for validation.
+        The model will set apart this fraction and will evaluate the loss and
+        any model metrics on this data at the end of each epoch.
+    **kwargs (dict): Other parameters to be passed to `lightgbm::lightgbm`
+        function.
+
+Returns:
+    SITSMachineLearningMethod: Model fitted to input data (to be passed to
+        `sits_classify`).
+
+Examples:
+    from pysits import *
+
+    # Example of training a model for time series classification
+    # Retrieve the samples for Mato Grosso
+    # train a lightgbm model
+    lgb_model = sits_train(samples_modis_ndvi, ml_method=sits_lightgbm)
+    # select the NDVI band
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    # classify the point
+    point_class = sits_classify(data=point_ndvi, ml_model=lgb_model)
+    plot(point_class)
diff --git a/pysits/docs/content/sits_lighttae.md b/pysits/docs/content/sits_lighttae.md
new file mode 100644
index 0000000..9cb56c8
--- /dev/null
+++ b/pysits/docs/content/sits_lighttae.md
@@ -0,0 +1,89 @@
+Train a model using Lightweight Temporal Self-Attention Encoder
+
+Implementation of Light Temporal Attention Encoder (L-TAE) for satellite image
+time series. This is a lightweight version of the temporal attention encoder
+proposed by Garnot et al. For the TAE, please see `sits_tae`.
+TAE is a simplified version of the well-known self-attention architecture which
+is used in large language models. Its modified self-attention scheme that uses
+the input embeddings as values. TAE defines a single master query for each
+sequence, computed from the temporal average of the queries. This master query
+is compared to the sequence of keys to produce a single attention mask used to
+weight the temporal mean of values into a single feature vector.
+The lightweight version of TAE further simplifies the TAE model. It defines
+master query of each head as a model parameter instead of the results of a
+linear layer, as is done it TAE. The authors argue that such simplification
+reduces the number of parameters, while the lack of flexibility is compensated
+by the larger number of available heads.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    samples_validation (SITSTimeSeriesModel): Time series with the validation
+        samples. If `samples_validation` parameter is provided,
+        `validation_split` is ignored.
+    epochs (int): Number of iterations to train the model (min = 1, max =
+        20000).
+    batch_size (int): Number of samples per gradient update (min = 16, max =
+        2048).
+    validation_split (float): Fraction of training data to be used as
+        validation data.
+    optimizer: Optimizer function to be used.
+    opt_hparams (dict): Hyperparameters for optimizer: `lr` : Learning rate of
+        the optimizer `eps`: Term added to the denominator to improve
+        numerical stability. `weight_decay`: L2 regularization rate.
+    lr_decay_epochs (int): Number of epochs to reduce learning rate.
+    lr_decay_rate (float): Decay factor for reducing learning rate.
+    patience (int): Number of epochs without improvements until training stops.
+    min_delta (float): Minimum improvement in loss function to reset the
+        patience counter.
+    seed (int): Seed for random values.
+    verbose (bool): Verbosity mode. Default is False.
+
+Returns:
+    R: A fitted model to be used for classification of data cubes.
+
+Notes:
+    `sits` provides a set of default values for all classification models.
+    These settings have been chosen based on testing by the authors.
+    Nevertheless, users can control all parameters for each model. Novice users
+    can rely on the default values, while experienced ones can fine-tune deep
+    learning models using `sits_tuning`.
+    This function is based on the paper by Vivien Garnot referenced below and
+    code available on github at https://github.com/VSainteuf/lightweight-
+    temporal-attention-pytorch If you use this method, please cite the original
+    TAE and the LTAE paper.
+    We also used the code made available by Maja Schneider in her work with
+    Marco K\u00f6rner referenced below and available at
+    https://github.com/maja601/RC2020-psetae.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a lightTAE model
+    torch_model = sits_train(samples_modis_ndvi, sits_lighttae())
+    # plot the model
+    plot(torch_model)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=torch_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_list_collections.md b/pysits/docs/content/sits_list_collections.md
new file mode 100644
index 0000000..5def432
--- /dev/null
+++ b/pysits/docs/content/sits_list_collections.md
@@ -0,0 +1,18 @@
+List the cloud collections supported by sits
+
+Prints the collections available in each cloud service supported by sits.
+Users can select to get information only for a single service by using the
+`source` parameter.
+
+Args:
+    source (str): Data source to be shown in detail.
+
+Returns:
+    None: Prints collections available in each cloud service supported by
+        sits.
+
+Examples:
+    from pysits import *
+
+    # show the names of the collections supported by SITS
+    sits_list_collections()
diff --git a/pysits/docs/content/sits_merge.md b/pysits/docs/content/sits_merge.md
new file mode 100644
index 0000000..1cd32c4
--- /dev/null
+++ b/pysits/docs/content/sits_merge.md
@@ -0,0 +1,44 @@
+Merge two data sets (time series or cubes)
+
+To merge two series, we consider that they contain different attributes but
+refer to the same data cube and spatiotemporal location. This function is
+useful for merging different bands of the same location. For example, one may
+want to put the raw and smoothed bands for the same set of locations in the
+same data set.
+In the case of data cubes, the function merges the images based on the
+following conditions:
+1. If the two cubes have different bands but compatible timelines, the bands
+   are combined, and the timeline is adjusted to overlap. To create the
+   overlap, we align the timelines like a "zipper": for each interval defined
+   by a pair of consecutive dates in the first timeline, we include matching
+   dates from the second timeline. If the second timeline has multiple dates in
+   the same interval, only the minimum date is kept. This ensures the final
+   timeline avoids duplicates and is consistent. This is useful when merging
+   data from different sensors (e.g., Sentinel-1 with Sentinel-2).
+2. If the bands are the same, the cube will have the combined timeline of both
+   cubes. This is useful for merging data from the same sensors from different
+   satellites (e.g., Sentinel-2A with Sentinel-2B).
+3. otherwise, the function will produce an error.
+
+Args:
+    data1 (SITSTimeSeriesModel | SITSCubeModel): Time series or data cube.
+    data2 (SITSTimeSeriesModel | SITSCubeModel): Time series or data cube.
+    suffix (list[str]): If data1 and data2 have duplicate bands, this
+        suffix will be added.
+    **kwargs (dict): Additional parameters.
+
+Returns:
+    SITSFrame: merged data sets (time series or data cube).
+
+Examples:
+    from pysits import *
+
+    # Retrieve a time series with values of NDVI
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    point_evi = sits_select(point_mt_6bands, bands="EVI")
+
+    # Merge time series
+    point_ndvi_evi = sits_merge(point_ndvi, point_evi)
+
+    # Plot the two points to see the smoothing effect
+    plot(point_ndvi_evi)
diff --git a/pysits/docs/content/sits_mixture_model.md b/pysits/docs/content/sits_mixture_model.md
new file mode 100644
index 0000000..4022a58
--- /dev/null
+++ b/pysits/docs/content/sits_mixture_model.md
@@ -0,0 +1,100 @@
+Multiple endmember spectral mixture analysis
+
+Create a multiple endmember spectral mixture analyses fractions images. We use
+the non-negative least squares (NNLS) solver to calculate the fractions of each
+endmember. The NNLS was implemented by Jakob Schwalb-Willmann in RStoolbox
+package (licensed as GPL>=3).
+
+Args:
+    data (SITSCubeModel | SITSTimeSeriesModel): A data cube or a set of
+        sample time series.
+    endmembers (pandas.DataFrame | str | pathlib.Path): Reference spectral
+        endmembers (see details below).
+    rmse_band (bool): Whether the error associated with the linear model
+        should be generated. If `True`, a new band with errors for each
+        pixel is generated using the root mean square measure (RMSE).
+        Default is `True`.
+    multicores (int): Number of cores to be used for generate the mixture
+        model.
+    progress (bool): Show progress bar? Default is `True`.
+    memsize (int): Memory available for the mixture model (in GB).
+    output_dir (str | pathlib.Path): Directory for output images.
+    **kwargs (dict): Parameters for specific functions.
+
+Returns:
+    SITSFrame: In case of a cube, a data cube with the fractions of each
+        endmember. The sum of all fractions is restricted to 1 (scaled from
+        0 to 10000), corresponding to the abundance of the endmembers in the
+        pixels. In case of a set of sample time series, the time series with
+        the values corresponding to each fraction.
+
+Notes:
+    Many pixels in images of medium-resolution satellites such as Landsat or
+    Sentinel-2 contain a mixture of spectral responses of different land cover
+    types. In many applications, it is desirable to obtain the proportion of a
+    given class inside a mixed pixel. For this purpose, the literature proposes
+    mixture models; these models represent pixel values as a combination of
+    multiple pure land cover types. Assuming that the spectral response of pure
+    land cover classes (called endmembers) is known, spectral mixture analysis
+    derives new bands containing the proportion of each endmember inside a
+    pixel.
+    The `endmembers` parameter should be a `pandas.DataFrame`, csv or a
+    shapefile. `endmembers` parameter must have the following columns: `type`,
+    which defines the endmembers that will be created and the columns
+    corresponding to the bands that will be used in the mixture model. The band
+    values must follow the product scale. For example, in the case of
+    sentinel-2 images the bands should be in the range 0 to 1. See the
+    `example` in this documentation for more details.
+
+Examples:
+    from pysits import *
+    import pandas as pd
+    import tempfile
+    import os
+
+    # Create a sentinel-2 cube
+    s2_cube = sits_cube(
+        source="AWS",
+        collection="SENTINEL-2-L2A",
+        tiles="20LKP",
+        bands=["B02", "B03", "B04", "B8A", "B11", "B12", "CLOUD"],
+        start_date="2019-06-13",
+        end_date="2019-06-30"
+    )
+    # create a directory to store the regularized file
+    reg_dir = os.path.join(tempfile.gettempdir(), "mix_model")
+    os.makedirs(reg_dir, exist_ok=True)
+    # Cube regularization for 16 days and 160 meters
+    reg_cube = sits_regularize(
+        cube=s2_cube,
+        period="P16D",
+        res=160,
+        roi={
+            "lon_min": -65.54870165,
+            "lat_min": -10.63479162,
+            "lon_max": -65.07629670,
+            "lat_max": -10.36046639
+        },
+        multicores=2,
+        output_dir=reg_dir
+    )
+
+    # Create the endmembers tibble
+    em = pd.DataFrame({
+        "class": ["forest", "land", "water"],
+        "B02": [0.02, 0.04, 0.07],
+        "B03": [0.0352, 0.065, 0.11],
+        "B04": [0.0189, 0.07, 0.14],
+        "B8A": [0.28, 0.36, 0.085],
+        "B11": [0.134, 0.35, 0.004],
+        "B12": [0.0546, 0.18, 0.0026]
+    })
+
+    # Generate the mixture model
+    mm = sits_mixture_model(
+        data=reg_cube,
+        endmembers=em,
+        memsize=4,
+        multicores=2,
+        output_dir=tempfile.gettempdir()
+    )
diff --git a/pysits/docs/content/sits_mlp.md b/pysits/docs/content/sits_mlp.md
new file mode 100644
index 0000000..f1ac6b3
--- /dev/null
+++ b/pysits/docs/content/sits_mlp.md
@@ -0,0 +1,85 @@
+Train multi-layer perceptron models using torch
+
+Use a multi-layer perceptron algorithm to classify data. This function uses
+the R "torch" and "luz" packages. Please refer to the documentation of those
+package for more details.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    samples_validation (SITSTimeSeriesModel): Time series with the
+        validation samples. If the `samples_validation` parameter is
+        provided, the `validation_split` parameter is ignored.
+    layers (list[int]): Number of hidden nodes in each layer.
+    dropout_rates (list[float]): Dropout rates (0,1) for each layer.
+    optimizer: Optimizer function to be used.
+    opt_hparams (dict): Hyperparameters for optimizer: lr : Learning rate
+        of the optimizer eps: Term added to the denominator to improve
+        numerical stability.. weight_decay: L2 regularization
+    epochs (int): Number of iterations to train the model.
+    batch_size (int): Number of samples per gradient update.
+    validation_split (float): Number between 0 and 1. Fraction of the
+        training data for validation. The model will set apart this
+        fraction and will evaluate the loss and any model metrics on this
+        data at the end of each epoch.
+    patience (int): Number of epochs without improvements until training
+        stops.
+    min_delta (float): Minimum improvement in loss function to reset the
+        patience counter.
+    seed (int): Seed for random values.
+    verbose (bool): Verbosity mode (True/False). Default is False.
+
+Returns:
+    SITSMachineLearningMethod: A torch mlp model to be used for
+        classification.
+
+Notes:
+    `sits` provides a set of default values for all classification models.
+    These settings have been chosen based on testing by the authors.
+    Nevertheless, users can control all parameters for each model. Novice
+    users can rely on the default values, while experienced ones can
+    fine-tune deep learning models using `sits_tuning`.
+    The default parameters for the MLP have been chosen based on the work by
+    Wang et al. 2017 that takes multilayer perceptrons as the baseline for
+    time series classifications: (a) Three layers with 512 neurons each,
+    specified by the parameter `layers`; (b) dropout rates of 10 (c) the
+    "optimizer_adam" as optimizer (default value); (d) a number of training
+    steps (`epochs`) of 100; (e) a `batch_size` of 64, which indicates how
+    many time series are used for input at a given steps; (f) a validation
+    percentage of 20 will be randomly set side for validation. (g) The
+    "relu" activation function.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create an MLP model
+    torch_model = sits_train(
+        samples_modis_ndvi,
+        sits_mlp(epochs=20, verbose=True)
+    )
+    # plot the model
+    plot(torch_model)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=torch_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_model_export.md b/pysits/docs/content/sits_model_export.md
new file mode 100644
index 0000000..1139189
--- /dev/null
+++ b/pysits/docs/content/sits_model_export.md
@@ -0,0 +1,20 @@
+Export classification models
+
+Given a trained machine learning or deep learning model, exports the
+model as an object for further exploration outside the `sits` package.
+
+Args:
+    ml_model (SITSMachineLearningMethod): A trained machine learning
+        model.
+
+Returns:
+    None: The model in the original format of the machine learning or
+        deep learning package.
+
+Examples:
+    from pysits import *
+
+    # create a classification model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # export the model
+    rfor_object = sits_model_export(rfor_model)
diff --git a/pysits/docs/content/sits_mosaic.md b/pysits/docs/content/sits_mosaic.md
new file mode 100644
index 0000000..5f4ca0b
--- /dev/null
+++ b/pysits/docs/content/sits_mosaic.md
@@ -0,0 +1,80 @@
+Mosaic classified cubes
+
+Creates a mosaic of all tiles of a data cube. Mosaics can be created from
+both regularized ARD images or from classified maps. In the case of ARD
+images, a mosaic will be produce for each band/date combination. It is
+better to first regularize the data cubes and then use `sits_mosaic`.
+
+Args:
+    cube (SITSCubeModel): A data cube.
+    crs (str | int): A target coordinate reference system of raster
+        mosaic. The provided crs could be a string (e.g, "EPSG:4326" or a
+        proj4string), or an EPSG code number (e.g. 4326). Default is
+        "EPSG:3857" - WGS 84 / Pseudo-Mercator.
+    roi (dict | str | pathlib.Path | geopandas.GeoDataFrame): Region of
+        interest (see below).
+    multicores (int): Number of cores that will be used to crop the images
+        in parallel.
+    output_dir (str | pathlib.Path): Directory for output images.
+    res (float): Spatial resolution of the mosaic. Default is None.
+    version (str): Version of resulting image (in the case of multiple
+        tests).
+    progress (bool): Show progress bar? Default is True.
+
+Returns:
+    SITSCubeModel: a data cube with only one tile.
+
+Notes:
+    To define a `roi` use one of:
+    - A path to a shapefile with polygons;
+    - A `geopandas.GeoDataFrame`;
+    - A named `dict` (`"lon_min"`, `"lat_min"`, `"lon_max"`,
+      `"lat_max"`) in WGS84;
+    - A named `dict` (`"xmin"`, `"xmax"`, `"ymin"`, `"ymax"`) with XY
+      coordinates.
+    The user should specify the CRS of the mosaic. We use "EPSG:3857"
+    (Pseudo-Mercator) as the default.
+
+Examples:
+    from pysits import *
+    import tempfile
+    from shapely.geometry import Polygon
+    import geopandas as gpd
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.mkdtemp()
+    )
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.mkdtemp())
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.mkdtemp()
+    )
+    # create roi
+    roi = gpd.GeoSeries(
+        [Polygon([
+            (-55.64768, -11.68649),
+            (-55.69654, -11.66455),
+            (-55.62973, -11.61519),
+            (-55.64768, -11.68649)
+        ])],
+        crs="EPSG:4326"
+    )
+    # crop and mosaic classified image
+    mosaic_cube = sits_mosaic(
+        cube=label_cube,
+        roi=roi,
+        crs="EPSG:4326",
+        output_dir=tempfile.mkdtemp()
+    )
diff --git a/pysits/docs/content/sits_parallel.md b/pysits/docs/content/sits_parallel.md
new file mode 100644
index 0000000..d3ef150
--- /dev/null
+++ b/pysits/docs/content/sits_parallel.md
@@ -0,0 +1,38 @@
+Configure or query sits parallel processing
+
+`sits_parallel` starts, restarts, stops, or gets a persistent `PSOCK`
+cluster used by `sits` functions that support parallel processing.
+To stop the cluster and free resources, call `sits_parallel(workers = 0)`
+(recommended). Calling with `workers = 1` has the same effect (parallel
+disabled).
+When called with no arguments, `sits_parallel()` returns the current
+cluster object. If no cluster is active, it returns `None`.
+
+Args:
+    workers (int): Number of workers to use. - `workers >= 2`: start or
+        restart a `PSOCK` cluster with `workers` workers. -
+        `workers <= 1`: stop any active cluster.
+    log (bool): If `True`, enables worker log/debug mode.
+    output_dir (str | pathlib.Path): Output directory where log files are
+        written when `log = True`.
+
+Returns:
+    None: If called with no arguments, returns the current `parallel`
+        cluster object or `None` if no cluster is active. If called with
+        `workers`, returns `None`.
+
+Notes:
+    This function is intended for long pipelines and production
+    environments where repeatedly creating and stopping clusters inside
+    each `sits` call is expensive. After `sits_parallel(workers = N)` is
+    called, `sits` functions can reuse the same cluster across multiple
+    calls.
+    When `workers >= 2`, worker processes inherit the current library
+    paths (`.libPaths()`) and selected environment variables required for
+    data access (for example, variables starting with `AWS_`).
+    When the streaming GPU pipeline is enabled
+    (`SITS_GPU_PIPELINE=stream`), `sits` functions attach the pipeline's
+    read and write stages to this cluster instead of creating a private
+    worker pool. Reads use at most `workers - 1` nodes and writes use one
+    node, so for a classification with `multicores` read slots start the
+    cluster with `workers = multicores + 1`.
diff --git a/pysits/docs/content/sits_patterns.md b/pysits/docs/content/sits_patterns.md
new file mode 100644
index 0000000..9634921
--- /dev/null
+++ b/pysits/docs/content/sits_patterns.md
@@ -0,0 +1,25 @@
+Find temporal patterns associated to a set of time series
+
+This function takes a set of time series samples as input estimates a set of
+patterns. The patterns are calculated using a GAM model. The idea is to use a
+formula of type y ~ s(x), where x is a temporal reference and y if the value of
+the signal. For each time, there will be as many predictions as there are
+sample values. The GAM model predicts a suitable approximation that fits the
+assumptions of the statistical model, based on a smooth function.
+This method is based on the "createPatterns" method of the R dtwSat package,
+which is also described in the reference paper.
+
+Args:
+    data (SITSTimeSeriesModel): Time series.
+    freq (int): Interval in days for estimates.
+    formula: Formula to be applied in the estimate.
+    **kwargs (dict): Any additional parameters.
+
+Returns:
+    SITSTimeSeriesPatternsModel: Time series with patterns.
+
+Examples:
+    from pysits import *
+
+    patterns = sits_patterns(cerrado_2classes)
+    plot(patterns)
diff --git a/pysits/docs/content/sits_pre_train.md b/pysits/docs/content/sits_pre_train.md
new file mode 100644
index 0000000..c6e3001
--- /dev/null
+++ b/pysits/docs/content/sits_pre_train.md
@@ -0,0 +1,51 @@
+Pre-train deep learning models for sits
+
+Runs self-supervised and pre-training for Earth observation time series
+using `sits`. The result of the function is an embedding model that can
+be used to produce embeddings from Earth observation data cubes.
+The function is a thin wrapper that validates inputs and dispatches to a
+user-selected deep learning pre-training method. The method is
+responsible for preparing training data, fitting the model, and
+returning an encoder that can be used later by `sits_encode`.
+Pre-training methods are created by factory functions available. The
+package offers three self-supervised learning methods:
+- `sits_ssl_lejepa` based on the LeJEPA method proposed by Balestriero
+  and LeCun (2025);
+- `sits_ssl_vicreg` based on the VICReg method proposed by Bardes,
+  Ponce, and LeCun (2022);
+- `sits_ssl_mae` based on masked autoencoders adapted to image time
+  series, similar to Tseng et al. (2024).
+Two supervised learning methods are available:
+- `sits_contrastive_learning` that adapts the supervised contrastive
+  learning proposed by Khosla et al.(2020) to image time series.
+- `sits_barlow_twins` that is a version for labelled time series of the
+  approach proposed by Zbontar et al.(2021)
+These factories return a function (closure) that implements the full
+pre-training procedure when called with `samples`.
+
+Args:
+    samples (SITSTimeSeriesModel): Time-series samples. Labels are
+        optional and may or may not be used depending on the selected
+        pre-training method.
+    rl_method (SITSRepresentationLearningMethod): A pre-training
+        representation learning method used to build an encoder that
+        generates embeddings (e.g., `sits_ssl_lejepa()` or
+        `sits_contrastive_learning()`). It must be a function that
+        takes `samples` and returns an encoder.
+
+Returns:
+    SITSRepresentationLearningMethod: A pre-trained deep learning
+    encoder and the metadata required for subsequent encoding (e.g.,
+    band order, feature naming, and normalization statistics, when
+    applicable).
+
+Examples:
+    from pysits import *
+
+    mae_model = sits_pre_train(
+        samples=samples_modis_ndvi,
+        rl_method=sits_ssl_mae(
+            encoder_model=sits_tempcnn(),
+            mask_ratio=0.5
+        )
+    )
diff --git a/pysits/docs/content/sits_pred_features.md b/pysits/docs/content/sits_pred_features.md
new file mode 100644
index 0000000..e563c84
--- /dev/null
+++ b/pysits/docs/content/sits_pred_features.md
@@ -0,0 +1,21 @@
+Obtain numerical values of predictors for time series samples
+
+Predictors are X-Y values required for machine learning algorithms,
+organized as a data table where each row corresponds to a training
+sample. The first two columns of the predictors table are categorical
+("label_id" and "label"). The other columns are the values of each band
+and time, organized first by band and then by time. This function
+returns the numeric values associated to each sample.
+
+Args:
+    pred (pandas.DataFrame): X-Y predictors, with one row per sample.
+
+Returns:
+    SITSFrame: The Y predictors for the sample, with one row per
+        sample.
+
+Examples:
+    from pysits import *
+
+    pred = sits_predictors(samples_modis_ndvi)
+    features = sits_pred_features(pred)
diff --git a/pysits/docs/content/sits_pred_normalize.md b/pysits/docs/content/sits_pred_normalize.md
new file mode 100644
index 0000000..49fcd21
--- /dev/null
+++ b/pysits/docs/content/sits_pred_normalize.md
@@ -0,0 +1,21 @@
+Normalize predictor values
+
+Most machine learning algorithms require data to be normalized. This
+applies to the "SVM" method and to all deep learning ones. To normalize
+the predictors, it is required that the statistics per band for each
+sample have been obtained by the `sits_stats` function.
+
+Args:
+    pred (pandas.DataFrame): X-Y predictors, with one row per sample.
+    stats (dict): Values of time series for Q02 and Q98 of the data
+        (two elements).
+
+Returns:
+    SITSFrame: Normalized predictor values.
+
+Examples:
+    from pysits import *
+
+    stats = sits_stats(samples_modis_ndvi)
+    pred = sits_predictors(samples_modis_ndvi)
+    pred_norm = sits_pred_normalize(pred, stats)
diff --git a/pysits/docs/content/sits_pred_references.md b/pysits/docs/content/sits_pred_references.md
new file mode 100644
index 0000000..ae6635d
--- /dev/null
+++ b/pysits/docs/content/sits_pred_references.md
@@ -0,0 +1,20 @@
+Obtain categorical id and predictor labels for time series samples
+
+Predictors are X-Y values required for machine learning algorithms,
+organized as a data table where each row corresponds to a training
+sample. The first two columns of the predictors table are categorical
+("label_id" and "label"). The other columns are the values of each band
+and time, organized first by band and then by time. This function
+returns the numeric values associated to each sample.
+
+Args:
+    pred (pandas.DataFrame): X-Y predictors with one row per sample.
+
+Returns:
+    list: Labels associated to training samples.
+
+Examples:
+    from pysits import *
+
+    pred = sits_predictors(samples_modis_ndvi)
+    ref = sits_pred_references(pred)
diff --git a/pysits/docs/content/sits_pred_sample.md b/pysits/docs/content/sits_pred_sample.md
new file mode 100644
index 0000000..08b11d7
--- /dev/null
+++ b/pysits/docs/content/sits_pred_sample.md
@@ -0,0 +1,20 @@
+Obtain a fraction of the predictors data frame
+
+Many machine learning algorithms (especially deep learning) use part of the
+original samples as test data to adjust its hyperparameters and to find an
+optimal point of convergence using gradient descent. This function extracts a
+fraction of the predictors to serve as test values for the deep learning
+algorithm.
+
+Args:
+    pred (pandas.DataFrame): X-Y predictors, one row per sample.
+    frac (float): Fraction of the X-Y predictors to be extracted.
+
+Returns:
+    SITSFrame: The chosen fraction of the X-Y predictors.
+
+Examples:
+    from pysits import *
+
+    pred = sits_predictors(samples_modis_ndvi)
+    pred_frac = sits_pred_sample(pred, frac=0.5)
diff --git a/pysits/docs/content/sits_predictors.md b/pysits/docs/content/sits_predictors.md
new file mode 100644
index 0000000..e8f634c
--- /dev/null
+++ b/pysits/docs/content/sits_predictors.md
@@ -0,0 +1,13 @@
+Obtain predictors for time series samples
+
+Predictors are X-Y values required for machine learning algorithms,
+organized as a data table where each row corresponds to a training
+sample. The first two columns of the predictors table are categorical
+(`label_id` and `label`). The other columns are the values of each band
+and time, organized first by band and then by time.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series samples.
+
+Returns:
+    SITSFrame: The predictors for the sample, with one row per sample.
diff --git a/pysits/docs/content/sits_reclassify.md b/pysits/docs/content/sits_reclassify.md
new file mode 100644
index 0000000..c07a943
--- /dev/null
+++ b/pysits/docs/content/sits_reclassify.md
@@ -0,0 +1,152 @@
+Reclassify a classified cube
+
+Reclassify a classification cube using a set of named expressions.
+For classified cubes, expressions relabel pixels based on logical
+conditions that may combine information from the classified cube and
+an optional mask cube.
+For probability cubes and probability vector cubes, expressions are
+used to group input labels into new labels by aggregating
+probabilities (summing probabilities of the selected input labels).
+
+Args:
+    cube (SITSCubeModel): Image cube to be reclassified (a classified
+        cube, probability cube, or probability vector cube).
+    mask (SITSCubeModel): Image cube with additional information to be
+        used in expressions. Used only for classified cube
+        reclassification.
+    rules (dict): Expressions to be evaluated. For classified cubes,
+        expressions must evaluate to a boolean and may refer to `cube`
+        and `mask`. For probability cubes and probability vector
+        cubes, each named rule selects one or more input labels (for
+        example using `cube %in% c(...)`). The probabilities of the
+        selected labels are summed to produce the new label given by
+        the rule name.
+    exclude_mask_na (bool): Should cube pixels be set to `None` when
+        `None` values are found in mask pixels? (default `True`).
+        Used only for classified cubes.
+    memsize (int): Memory available for processing in GB (min = 1,
+        max = 16384).
+    multicores (int): Number of cores to be used for processing
+        (min = 1, max = 2048).
+    output_dir (str | pathlib.Path): Directory where files will be
+        saved.
+    version (str): Version of resulting image.
+    progress (bool): Show progress bar?
+    **kwargs (dict): Other parameters for specific methods.
+
+Returns:
+    SITSCubeModel: An object of the same type as `cube`: a classified
+        cube for label cubes, or a probability cube and probability
+        vector cube for probability cubes and probability vector
+        cubes, respectively.
+
+Notes:
+    For classified cubes, reclassification changes the class assigned
+    to each pixel based on user-defined rules. Users should refer to
+    `cube` and `mask` to construct logical expressions. Expressions
+    are evaluated sequentially on the original classified values;
+    later rules override earlier ones.
+    For probability cubes and probability vector cubes,
+    reclassification is intended to group classes by combining
+    probabilities. Each named rule defines a new output label. For
+    each pixel, the probabilities of the selected input labels are
+    summed and assigned to the corresponding output label. Rules are
+    evaluated on the original probability layers.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # Example for probs_cube: group labels by summing probabilities
+
+    # Train a model
+    rf_model = sits_train(samples_modis_ndvi, ml_method=sits_rfor)
+
+    # Open a cube
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # Classify cube
+    probs_cube = sits_classify(
+        data=cube,
+        ml_model=rf_model,
+        output_dir=tempfile.gettempdir(),
+        version="classify"
+    )
+
+    # Reclassify probs_cube
+    probs_nat_veg = sits_reclassify(
+        cube=probs_cube,
+        rules={
+            "Cerrado": 'cube %in% c("Cerrado", "Forest")'
+        },
+        output_dir=tempfile.gettempdir()
+    )
+    plot(probs_nat_veg)
+
+    # Example for label cube: replacement of labels
+
+    # Open mask map
+    data_dir = r_package_dir("extdata/raster/prodes", package="sits")
+    prodes2021 = sits_cube(
+        source="USGS",
+        collection="LANDSAT-C2L2-SR",
+        data_dir=data_dir,
+        parse_info=[
+            "X1", "X2", "tile", "start_date", "end_date",
+            "band", "version"
+        ],
+        bands="class",
+        version="v20220606",
+        labels={
+            "1": "Forest", "2": "Water", "3": "NonForest",
+            "4": "NonForest2", "6": "d2007", "7": "d2008",
+            "8": "d2009", "9": "d2010", "10": "d2011",
+            "11": "d2012", "12": "d2013", "13": "d2014",
+            "14": "d2015", "15": "d2016", "16": "d2017",
+            "17": "d2018", "18": "r2010", "19": "r2011",
+            "20": "r2012", "21": "r2013", "22": "r2014",
+            "23": "r2015", "24": "r2016", "25": "r2017",
+            "26": "r2018", "27": "d2019", "28": "r2019",
+            "29": "d2020", "31": "r2020", "32": "Clouds2021",
+            "33": "d2021", "34": "r2021"
+        },
+        progress=False
+    )
+
+    # Open classification map
+    data_dir = r_package_dir("extdata/raster/classif", package="sits")
+    ro_class = sits_cube(
+        source="MPC",
+        collection="SENTINEL-2-L2A",
+        data_dir=data_dir,
+        parse_info=[
+            "X1", "X2", "tile", "start_date", "end_date",
+            "band", "version"
+        ],
+        bands="class",
+        labels={
+            "1": "ClearCut_Fire", "2": "ClearCut_Soil",
+            "3": "ClearCut_Veg", "4": "Forest"
+        },
+        progress=False
+    )
+
+    # Reclassify cube
+    ro_mask = sits_reclassify(
+        cube=ro_class,
+        mask=prodes2021,
+        rules={
+            "Old_Deforestation": 'mask %in% c("d2007", "d2008", "d2009", "d2010", "d2011", "d2012", "d2013", "d2014", "d2015", "d2016", "d2017", "d2018", "r2010", "r2011", "r2012", "r2013", "r2014", "r2015", "r2016", "r2017", "r2018", "d2019", "r2019", "d2020", "r2020", "r2021")',
+            "Water_Mask": 'mask == "Water"',
+            "NonForest_Mask": 'mask %in% c("NonForest", "NonForest2")'
+        },
+        memsize=4,
+        multicores=2,
+        output_dir=tempfile.gettempdir(),
+        version="ex_reclassify"
+    )
diff --git a/pysits/docs/content/sits_reduce.md b/pysits/docs/content/sits_reduce.md
new file mode 100644
index 0000000..f493340
--- /dev/null
+++ b/pysits/docs/content/sits_reduce.md
@@ -0,0 +1,63 @@
+Reduces a cube or samples from a summarization function
+
+Apply a temporal reduction from a named expression in a data cube or set
+of sample time series. In the case of data cubes, it materializes a new
+band in `output_dir`. The result will be a cube with only one date with
+the raster reduced from the function.
+
+`sits_reduce()` allows valid R expression to compute new bands. Use R
+syntax to pass an expression to this function. Besides arithmetic
+operators, you can use virtually any R function that can be applied to
+elements of a matrix. The provided functions must operate at line level
+in order to perform temporal reduction on a pixel.
+`sits_reduce()` Applies a function to each row of a matrix. In this
+matrix, each row represents a pixel and each column represents a single
+date. We provide some operations already implemented in the package to
+perform the reduce operation. See the list of available functions below:
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): set of sample time
+        series or data cube.
+    impute_fn: Imputation function to remove NA values.
+    memsize (int): Memory available for classification (in GB).
+    multicores (int): Number of cores to be used for classification.
+    output_dir (str | pathlib.Path): Directory where files will be
+        saved.
+    progress (bool): Show progress bar?
+    **kwargs (dict): Named expressions to be evaluated (see details).
+
+Returns:
+    SITSFrame: sample time series or data cube with new bands, produced
+        according to the requested expression.
+
+Notes:
+    The `t_sum()`, `t_std()`, `t_skewness()`, `t_kurtosis`, `t_mse`
+    indexes generate values greater than the limit of a two-byte
+    integer. Therefore, we save the images generated by these as
+    Float-32 with no scale.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # Reduce summarization function
+    point2 = sits_reduce(
+        sits_select(point_mt_6bands, bands="NDVI"),
+        NDVI_MEDIAN="t_median(NDVI)"
+    )
+
+    # Example of generation mean summarization from a cube
+    # Create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # Reduce NDVI band with mean function
+    cube_mean = sits_reduce(
+        data=cube,
+        NDVIMEAN="t_mean(NDVI)",
+        output_dir=tempfile.gettempdir()
+    )
diff --git a/pysits/docs/content/sits_reduce_imbalance.md b/pysits/docs/content/sits_reduce_imbalance.md
new file mode 100644
index 0000000..8efc09f
--- /dev/null
+++ b/pysits/docs/content/sits_reduce_imbalance.md
@@ -0,0 +1,51 @@
+Reduce imbalance in a set of samples
+
+Takes a set of samples with different labels and returns a new set. Deals
+with class imbalance using the synthetic minority oversampling technique
+(SMOTE) for oversampling. Undersampling is done using the SOM methods
+available in the `sits` package.
+
+Args:
+    samples (SITSTimeSeriesModel): Sample set to rebalance.
+    n_samples_over (int): Number of samples to oversample for classes with
+        samples less than this number.
+    n_samples_under (int): Number of samples to undersample for classes
+        with samples more than this number.
+    method (str): Method for oversampling (default = "smote").
+    multicores (int): Number of cores to process the data (default 2).
+
+Returns:
+    SITSTimeSeriesModel: A set of samples with reduced imbalance.
+
+Notes:
+    Many training samples for Earth observation data analysis are
+    imbalanced. This situation arises when the distribution of samples
+    associated with each label is uneven. Sample imbalance is an
+    undesirable property of a training set. Reducing sample imbalance
+    improves classification accuracy.
+    The function `sits_reduce_imbalance` increases the number of samples of
+    least frequent labels, and reduces the number of samples of most
+    frequent labels. To generate new samples, `sits` uses the SMOTE method
+    that estimates new samples by considering the cluster formed by the
+    nearest neighbors of each minority label.
+    To perform undersampling, `sits_reduce_imbalance`) builds a SOM map for
+    each majority label based on the required number of samples. Each
+    dimension of the SOM is set to ceiling(sqrt(new_number_samples/4)) to
+    allow a reasonable number of neurons to group similar samples. After
+    calculating the SOM map, the algorithm extracts four samples per neuron
+    to generate a reduced set of samples that approximates the variation of
+    the original one. See also `sits_som_map`.
+
+Examples:
+    from pysits import *
+
+    # print the labels summary for a sample set
+    summary(samples_modis_ndvi)
+    # reduce the sample imbalance
+    new_samples = sits_reduce_imbalance(samples_modis_ndvi,
+        n_samples_over=200,
+        n_samples_under=200,
+        multicores=1
+    )
+    # print the labels summary for the rebalanced set
+    summary(new_samples)
diff --git a/pysits/docs/content/sits_regularize.md b/pysits/docs/content/sits_regularize.md
new file mode 100644
index 0000000..cb85b94
--- /dev/null
+++ b/pysits/docs/content/sits_regularize.md
@@ -0,0 +1,140 @@
+Build a regular data cube from an irregular one
+
+Produces regular data cubes for analysis-ready data (ARD) image collections.
+Analysis-ready data (ARD) collections available in AWS, MPC, USGS and DEAfrica
+are not regular in space and time. Bands may have different resolutions, images
+may not cover the entire time, and time intervals are not regular. For this
+reason, subsets of these collection need to be converted to regular data cubes
+before further processing and data analysis. This function requires users to
+include the cloud band in their ARD-based data cubes. This function uses the
+`gdalcubes` package.
+
+Args:
+    cube (SITSCubeModel): data cube whose observation period and/or spatial
+        resolution is not constant.
+    period (str): ISO8601-compliant time period for regular data cubes, with
+        number and unit, where "D", "M" and "Y" stand for days, month and
+        year; e.g., "P16D" for 16 days.
+    res (float): Spatial resolution of regularized images (in meters).
+    output_dir (str | pathlib.Path): Valid directory for storing regularized
+        images.
+    timeline (list[str]): User-defined timeline for regularized cube.
+    roi (dict | geopandas.GeoDataFrame | str | pathlib.Path): Region of
+        interest (see notes below).
+    crs (str | int): Coordinate Reference System (CRS) of the roi. (see
+        details below).
+    tiles (list[str]): Tiles to be produced.
+    grid_system (str): Grid system to be used for the output images.
+    multicores (int): Number of cores used for regularization; used for
+        parallel processing of input.
+    progress (bool): show progress bar?
+    **kwargs (dict): Additional parameters.
+
+Returns:
+    SITSCubeModel: A data cube with aggregated images.
+
+Notes:
+    The main `sits` classification workflow has the following steps:
+    1. `sits_cube`: selects a ARD image collection from a cloud provider.
+    2. `sits_cube_copy`: copies an ARD image collection from a cloud provider
+       to a local directory for faster processing.
+    3. `sits_regularize`: create a regular data cube from an ARD image
+       collection.
+    4. `sits_apply`: create new indices by combining bands of a regular data
+       cube (optional).
+    5. `sits_get_data`: extract time series from a regular data cube based on
+       user-provided labelled samples.
+    6. `sits_train`: train a machine learning model based on image time series.
+    7. `sits_classify`: classify a data cube using a machine learning model and
+       obtain a probability cube.
+    8. `sits_smooth`: post-process a probability cube using a spatial smoother
+       to remove outliers and increase spatial consistency.
+    9. `sits_label_classification`: produce a classified map by selecting the
+       label with the highest probability from a smoothed cube.
+    The regularization operation converts subsets of image collections
+    available in cloud providers into regular data cubes. It is an essential
+    part of the `sits` workflow. The input to `sits_regularize` should be an
+    ARD cube which includes the cloud band. The aggregation method used in
+    `sits_regularize` sorts the images based on cloud cover, putting images
+    with the least clouds at the top of the stack. Once the stack of images is
+    sorted, the method uses the first valid value to create the temporal
+    aggregation.
+    The "period" parameter is mandatory, and defines the time interval between
+    two images of the regularized cube. When combining Sentinel-1A and
+    Sentinel-1B images, experiments show that a 16-day period ("P16D") are a
+    good default. Landsat images require a longer period of one to three
+    months.
+    By default, the date of the first image of the input cube is taken as the
+    starting date for the regular cube. In many situations, users may want to
+    pre-define the required times using the "timeline" parameter. The
+    "timeline" parameter, if used, must contain a set of dates which are
+    compatible with the input cube.
+    To define a `roi` use one of:
+    - A path to a shapefile with polygons;
+    - A `geopandas.GeoDataFrame`;
+    - A `SpatExtent` object from `terra` package;
+    - A `dict` (`"lon_min"`, `"lat_min"`, `"lon_max"`, `"lat_max"`) in
+      WGS84;
+    - A `dict` (`"xmin"`, `"xmax"`, `"ymin"`, `"ymax"`) with XY
+      coordinates.
+    Defining a region of interest using `SpatExtent` or XY values not in WGS84
+    requires the `crs` parameter to be specified. `sits_regularize()` function
+    will crop the images that contain the region of interest(). NOTE: Make sure
+    to inform `roi` with valid geometries. `sits` will drop the use of `roi` if
+    it contains invalid geometries.
+    The optional `tiles` parameter indicates which tiles of the input cube will
+    be used for regularization. When `grid_system` is informed, `tiles` may be
+    combined with `roi` to further restrict which target grid tiles are
+    produced (only tiles that both intersect `roi` and are listed in `tiles`
+    are kept).
+    The `grid_system` parameter allows the user to reproject the files to a
+    grid system which is different from that used in the ARD image collection
+    of the could provider. Currently, the package supports the use of MGRS grid
+    system and those used by the Brazil Data Cube ("BDC_LG_V2" "BDC_MD_V2"
+    "BDC_SM_V2").
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # define a non-regular Sentinel-2 cube in AWS
+    s2_cube_open = sits_cube(
+        source="AWS",
+        collection="SENTINEL-2-L2A",
+        tiles=["20LKP", "20LLP"],
+        bands=["B8A", "CLOUD"],
+        start_date="2018-10-01",
+        end_date="2018-11-01"
+    )
+    # regularize the cube
+    rg_cube = sits_regularize(
+        cube=s2_cube_open,
+        period="P16D",
+        res=60,
+        multicores=2,
+        output_dir=tempfile.gettempdir()
+    )
+
+    ## Sentinel-1 SAR
+    roi = {
+        "lon_min": -50.410, "lon_max": -50.379,
+        "lat_min": -10.1910, "lat_max": -10.1573
+    }
+    s1_cube_open = sits_cube(
+        source="MPC",
+        collection="SENTINEL-1-GRD",
+        bands=["VV", "VH"],
+        orbit="descending",
+        roi=roi,
+        start_date="2020-06-01",
+        end_date="2020-09-28"
+    )
+    # regularize the cube
+    rg_cube = sits_regularize(
+        cube=s1_cube_open,
+        period="P12D",
+        res=60,
+        roi=roi,
+        multicores=2,
+        output_dir=tempfile.gettempdir()
+    )
diff --git a/pysits/docs/content/sits_resnet.md b/pysits/docs/content/sits_resnet.md
new file mode 100644
index 0000000..cb1527e
--- /dev/null
+++ b/pysits/docs/content/sits_resnet.md
@@ -0,0 +1,79 @@
+Train ResNet classification models
+
+Use a ResNet architecture for classifying image time series. The ResNet (or
+deep residual network) was proposed by a team in Microsoft Research for 2D
+image classification. ResNet tries to address the degradation of accuracy in a
+deep network. The idea is to replace a deep network with a combination of
+shallow ones. In the paper by Fawaz et al. (2019), ResNet was considered the
+best method for time series classification, using the UCR dataset. Please refer
+to the paper for more details.
+The R-torch version is based on the code made available by Zhiguang Wang,
+author of the original paper. The code was developed in python using keras.
+https://github.com/cauchyturing (repo:
+UCR_Time_Series_Classification_Deep_Learning_Baseline)
+The R-torch version also considered the code by Ignacio Oguiza, whose
+implementation is available at
+https://github.com/timeseriesAI/tsai/blob/main/tsai/models/ResNet.py.
+There are differences between Wang's Keras code and Oguiza torch code. In this
+case, we have used Wang's keras code as the main reference.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    samples_validation (SITSTimeSeriesModel): Time series with the
+        validation samples. If provided, `validation_split` is ignored.
+    blocks (list[int]): Number of 1D convolutional filters for each block
+        of three layers.
+    kernels (list[int]): Size of the 1D convolutional kernels.
+    epochs (int): Number of iterations to train the model, for each layer
+        of each block.
+    batch_size (int): Number of samples per gradient update.
+    validation_split (float): Fraction of training data to be used as
+        validation data.
+    optimizer: Optimizer function to be used.
+    opt_hparams (dict): Hyperparameters for optimizer: lr : Learning rate
+        of the optimizer eps: Term added to the denominator to improve
+        numerical stability. weight_decay: L2 regularization
+    lr_decay_epochs (int): Number of epochs to reduce learning rate.
+    lr_decay_rate (float): Decay factor for reducing learning rate.
+    patience (int): Number of epochs without improvements until training
+        stops.
+    min_delta (float): Minimum improvement in loss function to reset the
+        patience counter.
+    seed (int): Seed for random values.
+    verbose (bool): Verbosity mode. Default is `False`.
+
+Returns:
+    SITSMachineLearningMethod: A fitted model to be used for classification.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a ResNet model
+    torch_model = sits_train(samples_modis_ndvi, sits_resnet())
+    # plot the model
+    plot(torch_model)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=torch_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_rfor.md b/pysits/docs/content/sits_rfor.md
new file mode 100644
index 0000000..14b08e3
--- /dev/null
+++ b/pysits/docs/content/sits_rfor.md
@@ -0,0 +1,37 @@
+Train random forest models
+
+Use Random Forest algorithm to classify samples. This function is a
+front-end to the `randomForest` package. Please refer to the
+documentation in that package for more details.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    num_trees (int): Number of trees to grow. This should not be set to
+        too small a number, to ensure that every input row gets predicted
+        at least a few times (default: 100, min = 20).
+    mtry (int): Number of variables randomly sampled as candidates at each
+        split (default: None - use default value of
+        `randomForest::randomForest()` function, i.e.
+        `floor(sqrt(features))`).
+    classwt (dict): Assigns priors to classes, influencing the Gini index
+        for splitting. Note that this parameter affects the tree-building
+        process rather than just post-hoc voting.
+    **kwargs (dict): Other parameters to be passed to
+        `randomForest::randomForest` function.
+
+Returns:
+    SITSMachineLearningMethod: Model fitted to input data (to be passed to
+        `sits_classify`).
+
+Examples:
+    from pysits import *
+
+    # Example of training a model for time series classification
+    # Retrieve the samples for Mato Grosso
+    # train a random forest model
+    rf_model = sits_train(samples_modis_ndvi, ml_method=sits_rfor)
+    # classify the point
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    # classify the point
+    point_class = sits_classify(data=point_ndvi, ml_model=rf_model)
+    plot(point_class)
diff --git a/pysits/docs/content/sits_roi_to_tiles.md b/pysits/docs/content/sits_roi_to_tiles.md
new file mode 100644
index 0000000..6450552
--- /dev/null
+++ b/pysits/docs/content/sits_roi_to_tiles.md
@@ -0,0 +1,45 @@
+Find tiles of a given ROI and Grid System
+
+Given an ROI and grid system, this function finds the intersected tiles and
+returns them as a `geopandas.GeoDataFrame`.
+
+Args:
+    roi (dict | str | pathlib.Path | geopandas.GeoDataFrame): Region of
+        interest (see notes below).
+    crs (str): Coordinate Reference System (CRS) of the roi. (see details
+        below).
+    grid_system (str): Grid system to be used for the output images.
+        (Default is "MGRS")
+
+Returns:
+    SITSFrameSF: A `geopandas.GeoDataFrame` with the intersect tiles with
+        three columns tile_id, epsg, and the percentage of coverage area.
+
+Notes:
+    To define a `roi` use one of:
+    - A path to a shapefile with polygons;
+    - A `geopandas.GeoDataFrame`;
+    - A named `dict` (`"lon_min"`, `"lat_min"`, `"lon_max"`, `"lat_max"`) in
+      WGS84;
+    - A named `dict` (`"xmin"`, `"xmax"`, `"ymin"`, `"ymax"`) with XY
+      coordinates.
+    Defining a region of interest using XY values not in WGS84
+    requires the `crs` parameter to be specified.
+    The `grid_system` parameter allows the user to reproject the files to a
+    grid system which is different from that used in the ARD image collection
+    of the could provider. Currently, the package supports the use of MGRS grid
+    system and those used by the Brazil Data Cube ("BDC_LG_V2" "BDC_MD_V2"
+    "BDC_SM_V2").
+
+Examples:
+    from pysits import *
+
+    # Defining a ROI
+    roi = {
+        "lon_min": -64.037,
+        "lat_min": -9.644,
+        "lon_max": -63.886,
+        "lat_max": -9.389,
+    }
+    # Finding tiles
+    tiles = sits_roi_to_tiles(roi, grid_system="MGRS")
diff --git a/pysits/docs/content/sits_sample.md b/pysits/docs/content/sits_sample.md
new file mode 100644
index 0000000..f6b82db
--- /dev/null
+++ b/pysits/docs/content/sits_sample.md
@@ -0,0 +1,27 @@
+Sample a time series
+
+Takes samples from a set of time series and returns a new
+`SITSTimeSeriesModel`. For a given field as a group criterion, the
+result contains a percentage of the total number of samples per group.
+If frac > 1, all sampling will be done with replacement.
+
+Args:
+    data (SITSTimeSeriesModel): Set of time series.
+    frac (float): Percentage of samples to extract (range: 0.0 to 2.0,
+        default = 0.2).
+    oversample (bool): Oversample classes with small number of samples?
+
+Returns:
+    SITSTimeSeriesModel: A set of time series with a fixed quantity of
+        samples.
+
+Examples:
+    from pysits import *
+
+    # Retrieve a set of time series with 2 classes (cerrado_2classes)
+    # Print the labels of the resulting tibble
+    print(sits_labels_summary(cerrado_2classes))
+    # Sample by fraction
+    data_02 = sits_sample(cerrado_2classes, frac=0.2)
+    # Print the labels
+    print(sits_labels_summary(data_02))
diff --git a/pysits/docs/content/sits_sampling_design.md b/pysits/docs/content/sits_sampling_design.md
new file mode 100644
index 0000000..b9de38c
--- /dev/null
+++ b/pysits/docs/content/sits_sampling_design.md
@@ -0,0 +1,47 @@
+Allocation of sample size to strata
+
+Takes a class cube with different labels and allocates a number of sample
+sizes per strata to obtain suitable values of error-adjusted area,
+providing five allocation strategies.
+
+Args:
+    cube (SITSCubeModel): Classified data cube.
+    expected_ua (dict): Expected values of user's accuracy.
+    alloc_options (list[float]): Fixed sample allocation for rare
+        classes.
+    std_err (float): Standard error we would like to achieve.
+    rare_class_prop (float): Proportional area limit for rare classes.
+
+Returns:
+    SITSMatrix: options to decide allocation of sample size to each
+        class. This uses the same format as Table 5 of
+        Olofsson et al.(2014).
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # label the probability cube
+    label_cube = sits_label_classification(
+        probs_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # estimated UA for classes
+    expected_ua = {
+        "Cerrado": 0.75, "Forest": 0.9,
+        "Pasture": 0.8, "Soy_Corn": 0.8
+    }
+    sampling_design = sits_sampling_design(label_cube, expected_ua)
diff --git a/pysits/docs/content/sits_segment.md b/pysits/docs/content/sits_segment.md
new file mode 100644
index 0000000..d3f7e2b
--- /dev/null
+++ b/pysits/docs/content/sits_segment.md
@@ -0,0 +1,91 @@
+Segment an image
+
+Apply a spatial-temporal segmentation on a data cube based on a user defined
+segmentation function. The function applies the segmentation algorithm "seg_fn"
+to each tile. The output is a vector data cube, which is a data cube with an
+additional vector file in "geopackage" format.
+
+Args:
+    cube (SITSCubeModel): Regular data cube.
+    seg_fn: Function to apply the segmentation.
+    roi (dict | str | pathlib.Path | geopandas.GeoDataFrame): Region of
+        interest (see below).
+    impute_fn: Imputation function to remove NA values.
+    start_date (str): Start date for the segmentation.
+    end_date (str): End date for the segmentation.
+    memsize (int): Memory available for classification (in GB).
+    multicores (int): Number of cores to be used for classification.
+    output_dir (str | pathlib.Path): Directory for output file.
+    version (str): Version of the output (for multiple segmentations).
+    progress (bool): Show progress bar?
+    **kwargs (dict): Additional arguments.
+
+Returns:
+    SITSCubeModel: The segmentation as a vector data cube.
+
+Notes:
+    Segmentation requires the following steps:
+    1. Create a regular data cube with `sits_cube` and `sits_regularize`;
+    2. Run `sits_segment` to obtain a vector data cube with polygons that
+       define the boundary of the segments;
+    3. Classify the time series associated to the segments with
+       `sits_classify`, to get obtain a vector probability cube;
+    4. Use `sits_label_classification` to label the vector probability cube;
+    5. Display the results with `plot` or `sits_view`.
+    The "roi" parameter defines a region of interest. It can be a
+    `geopandas.GeoDataFrame`, a shapefile, or a bounding box `dict` with named
+    XY values ("xmin", "xmax", "ymin", "ymax") or named lat/long values
+    ("lon_min", "lat_min", "lon_max", "lat_max").
+    As of version 1.5.4, two segmentation functions are available. The
+    preferred option is `sits_snic`, which implements the Simple Non-Iterative
+    Clustering (SNIC) algorithm to generate compact and homogeneous superpixels
+    directly from uniformly distributed seeds. SNIC avoids the iterative
+    refinement step used in SLIC and is generally faster and more memory-
+    efficient, making it suitable for large multispectral or multitemporal data
+    cubes.
+    The previous function `sits_slic`, based on the Simple Linear Iterative
+    Clustering (SLIC) algorithm as adapted by Nowosad and Stepinski for
+    multispectral and multitemporal imagery, remains available but is now
+    deprecated and will be removed in a future release. SLIC clusters pixels
+    using spectral similarity and spatial–temporal proximity to produce nearly
+    uniform superpixels, but its iterative nature makes it less efficient for
+    large-scale Earth observation workflows.
+    The result of `sits_segment` is a data cube with an additional vector file
+    in the `geopackage` format. The location of the vector file is included in
+    the data cube in a new column, called `vector_info`.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    # create a data cube
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # segment the vector cube
+    segments = sits_segment(
+        cube=cube,
+        seg_fn=sits_snic(
+            grid_seeding="diamond",
+            spacing=15,
+            compactness=0.5,
+            padding=2
+        ),
+        output_dir=tempfile.mkdtemp()
+    )
+    # create a classification model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # classify the segments
+    seg_probs = sits_classify(
+        data=segments,
+        ml_model=rfor_model,
+        output_dir=tempfile.mkdtemp()
+    )
+    # label the probability segments
+    seg_label = sits_label_classification(
+        cube=seg_probs,
+        output_dir=tempfile.mkdtemp()
+    )
diff --git a/pysits/docs/content/sits_select.md b/pysits/docs/content/sits_select.md
new file mode 100644
index 0000000..cc6641f
--- /dev/null
+++ b/pysits/docs/content/sits_select.md
@@ -0,0 +1,38 @@
+Filter a data set for bands, tiles, and dates
+
+Filter the bands, tiles, dates and labels from a set of time series or
+from a data cube.
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): time series or data
+        cube.
+    bands (list[str]): names of the bands.
+    start_date (str): date in YYYY-MM-DD format: start date to be
+        filtered.
+    end_date (str): date in YYYY-MM-DD format: end date to be
+        filtered.
+    dates (list[str]): sparse dates to be selected.
+    labels (list[str]): sparse labels to be selected (only applied
+        for time series data).
+    tiles (list[str]): names of the tiles.
+    **kwargs (dict): additional parameters to be provided.
+
+Returns:
+    SITSFrame: time series or data cube.
+
+Examples:
+    from pysits import *
+
+    # Retrieve a set of time series with 2 classes
+    # (cerrado_2classes is available directly)
+    # Print the original bands
+    print(sits_bands(cerrado_2classes))
+    # Select only the NDVI band
+    data = sits_select(cerrado_2classes, bands=["NDVI"])
+    # Print the labels of the resulting tibble
+    print(sits_bands(data))
+    # select start and end date
+    point_2010 = sits_select(point_mt_6bands,
+        start_date="2000-09-13",
+        end_date="2017-08-29"
+    )
diff --git a/pysits/docs/content/sits_show_prediction.md b/pysits/docs/content/sits_show_prediction.md
new file mode 100644
index 0000000..24b3636
--- /dev/null
+++ b/pysits/docs/content/sits_show_prediction.md
@@ -0,0 +1,23 @@
+Shows the predicted labels for a classified time series
+
+This function takes a classified time series produced by a machine
+learning method and displays the result.
+
+Args:
+    class (SITSTimeSeriesModel): A time series that has been classified.
+
+Returns:
+    SITSFrame: Table with the columns "from", "to", "class".
+
+Examples:
+    from pysits import *
+
+    # Retrieve the samples for Mato Grosso
+    # train an SVM model
+    ml_model = sits_train(samples_modis_ndvi, ml_method=sits_svm)
+    # classify the point
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    point_class = sits_classify(
+        data=point_ndvi, ml_model=ml_model
+    )
+    sits_show_prediction(point_class)
diff --git a/pysits/docs/content/sits_slic.md b/pysits/docs/content/sits_slic.md
new file mode 100644
index 0000000..5397b23
--- /dev/null
+++ b/pysits/docs/content/sits_slic.md
@@ -0,0 +1,71 @@
+Segment an image using SLIC
+
+Apply a segmentation on a data cube using either the `supercells` or `snic`
+packages, depending on the chosen algorithm. As of version 1.5.4, two
+segmentation methods are supported. The recommended option is SNIC, implemented
+via the `snic` package, which applies a non-iterative clustering strategy to
+generate compact, homogeneous superpixels from uniformly distributed seeds
+(Achanta and Susstrunk, 2017). The alternative method uses the SLIC algorithm
+implemented in the `supercells` package, adapted for remote sensing data
+following Achanta et al. (2012). This SLIC variant is deprecated and will be
+removed in a future release. See references for more details.
+
+Args:
+    data (SITSMatrix): time series values.
+    step (int): distance (in number of cells) between initial supercells'
+        centers.
+    compactness (float): compactness value. Larger values cause clusters to
+        be more compact/even (square).
+    dist_fun (str): distance function. Currently implemented: `euclidean,
+        jsd, dtw`, and any distance function from the `philentropy` package.
+        See `philentropy::getDistMethods()`.
+    avg_fun (str): averaging function to calculate the values of the
+        supercells' centers. Accepts any fitting function (e.g., mean or
+        median) or one of internally implemented "mean" and "median".
+        Default: "median".
+    iter (int): number of iterations to create the output.
+    minarea (int): minimal size of a supercell (in cells).
+    verbose (bool): show the progress bar?
+
+Returns:
+    R: set of segments for a single tile.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    # create a data cube
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # segment the vector cube
+    segments = sits_segment(
+        cube=cube,
+        seg_fn=sits_snic(
+            grid_seeding="rectangular",
+            spacing=10,
+            compactness=0.3,
+            padding=0
+        ),
+        output_dir=tempfile.gettempdir(),
+        version="snic-demo"
+    )
+    # create a classification model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # classify the segments
+    seg_probs = sits_classify(
+        data=segments,
+        ml_model=rfor_model,
+        output_dir=tempfile.gettempdir(),
+        version="snic-demo"
+    )
+    # label the probability segments
+    seg_label = sits_label_classification(
+        cube=seg_probs,
+        output_dir=tempfile.gettempdir(),
+        version="snic-demo"
+    )
+    plot(seg_label)
diff --git a/pysits/docs/content/sits_smooth.md b/pysits/docs/content/sits_smooth.md
new file mode 100644
index 0000000..a9a756f
--- /dev/null
+++ b/pysits/docs/content/sits_smooth.md
@@ -0,0 +1,90 @@
+Smooth probability cubes with spatial predictors
+
+Takes a set of classified raster layers with probabilities and applies a
+Bayesian smoothing function.
+
+Args:
+    cube (SITSCubeModel): Probability data cube.
+    window_size (int): Size of the neighborhood (min = 3, max = 21).
+    neigh_fraction (float): Fraction of neighbors with high probabilities
+        to be used in Bayesian inference (min = 0.1, max = 1.0).
+    smoothness (int | list[int]): Estimated variance of logit of class
+        probabilities (Bayesian smoothing parameter) (min = 1, max = 200).
+    exclusion_mask (geopandas.GeoDataFrame | str | pathlib.Path): Areas to
+        be excluded from the classification process.
+    memsize (int): Memory available for classification in GB (min = 1,
+        max = 16384).
+    multicores (int): Number of cores to be used for classification
+        (min = 1, max = 2048).
+    output_dir (str | pathlib.Path): Valid directory for output file.
+    version (str): Version of the output.
+    progress (bool): Check progress bar?
+    **kwargs (dict): Other parameters for specific functions.
+
+Returns:
+    SITSCubeModel: A data cube.
+
+Notes:
+    The main `sits` classification workflow has the following steps:
+    1. `sits_cube`: selects a ARD image collection from a cloud provider.
+    2. `sits_cube_copy`: copies an ARD image collection from a cloud provider
+       to a local directory for faster processing.
+    3. `sits_regularize`: create a regular data cube from an ARD image
+       collection.
+    4. `sits_apply`: create new indices by combining bands of a regular data
+       cube (optional).
+    5. `sits_get_data`: extract time series from a regular data cube based on
+       user-provided labelled samples.
+    6. `sits_train`: train a machine learning model based on image time series.
+    7. `sits_classify`: classify a data cube using a machine learning model and
+       obtain a probability cube.
+    8. `sits_smooth`: post-process a probability cube using a spatial smoother
+       to remove outliers and increase spatial consistency.
+    9. `sits_label_classification`: produce a classified map by selecting the
+       label with the highest probability from a smoothed cube.
+    Machine learning algorithms rely on training samples that are derived from
+    \--picked by users to represent the desired output
+    classes. Given the presence of mixed pixels in images regardless of
+    resolution, and the considerable data variability within each class, these
+    classifiers often produce results with misclassified pixels.
+    Post-processing the results of `sits_classify` using `sits_smooth` reduces
+    salt-and-pepper and border effects. By minimizing noise, `sits_smooth`
+    brings a significant gain in the overall accuracy and interpretability of
+    the final output.
+    When applied to a `probs_cube`, `sits_smooth` uses a spatial window defined
+    by the `window_size` parameter to identify neighboring pixels. When applied
+    to a `probs_vector_cube`, the function uses segments to define neighbors.
+    All pixels within a segment are considered neighbors. Together with
+    `neigh_fraction`, this determines the fraction of pixels overlapped by the
+    segment that are used in the smoothing process.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create an xgboost model
+    xgb_model = sits_train(samples_modis_ndvi, sits_xgboost())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=xgb_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_snic.md b/pysits/docs/content/sits_snic.md
new file mode 100644
index 0000000..cc12fc5
--- /dev/null
+++ b/pysits/docs/content/sits_snic.md
@@ -0,0 +1,60 @@
+Segment an image using SNIC
+
+Apply a segmentation on a data cube based on the `snic` package. This is
+an adaptation and extension to remote sensing data of the SNIC
+superpixels algorithm proposed by Achanta and S\u00fcsstrunk (2017). See
+reference for more details.
+
+Args:
+    data (pandas.DataFrame): Time series.
+    grid_seeding (str): Method for grid seeding (one of "rectangular",
+        "diamond", "hexagonal", "random").
+    spacing (int): Distance (in number of cells) between initial
+        supercells' centers.
+    compactness (float): A compactness value. Larger values cause
+        clusters to be more compact/even (square).
+    padding (int): Distance (in pixels) from the image borders within
+        which no seeds are placed.
+
+Returns:
+    R: The segmentation function to be applied to a data cube.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    # create a data cube
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # segment the vector cube
+    segments = sits_segment(
+        cube=cube,
+        seg_fn=sits_snic(
+            grid_seeding="rectangular",
+            spacing=10,
+            compactness=0.5,
+            padding=5
+        ),
+        output_dir=tempfile.mkdtemp(),
+        version="snic-demo"
+    )
+    # create a classification model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # classify the segments
+    seg_probs = sits_classify(
+        data=segments,
+        ml_model=rfor_model,
+        output_dir=tempfile.mkdtemp(),
+        version="snic-demo"
+    )
+    # label the probability segments
+    seg_label = sits_label_classification(
+        cube=seg_probs,
+        output_dir=tempfile.mkdtemp(),
+        version="snic-demo"
+    )
+    plot(seg_label)
diff --git a/pysits/docs/content/sits_som_clean_samples.md b/pysits/docs/content/sits_som_clean_samples.md
new file mode 100644
index 0000000..701a30a
--- /dev/null
+++ b/pysits/docs/content/sits_som_clean_samples.md
@@ -0,0 +1,45 @@
+Cleans the samples based on SOM map information
+
+`sits_som_clean_samples()` evaluates the quality of the samples based on the
+results of the SOM map.
+
+Args:
+    som_map: returned by `sits_som_map`.
+    prior_threshold (float): threshold of conditional probability (frequency
+        of samples assigned to the same SOM neuron).
+    posterior_threshold (float): threshold of posterior probability
+        (influenced by the SOM neighborhood).
+    keep (list[str]): which types of evaluation to be maintained in the data.
+
+Returns:
+    SITSTimeSeriesModel: samples with two additional columns. The first
+        indicates if each sample is clean, should be analyzed or should be
+        removed. The second is the posterior probability of the sample. The
+        "keep" parameter indicates which.
+
+Notes:
+    The algorithm identifies noisy samples, using `prior_threshold` for the
+    prior probability and `posterior_threshold` for the posterior probability.
+    Each sample receives an evaluation tag, according to the following rule:
+    (a) If the prior probability is < `prior_threshold`, the sample is tagged
+    as "remove"; (b) If the prior probability is >= `prior_threshold` and the
+    posterior probability is >=`posterior_threshold`, the sample is tagged as
+    "clean"; (c) If the prior probability is >= `posterior_threshold` and the
+    posterior probability is < `posterior_threshold`, the sample is tagged as
+    "analyze" for further inspection. The user can define which tagged samples
+    will be returned using the "keep" parameter, with the following options:
+    "clean", "analyze", "remove".
+
+Examples:
+    from pysits import *
+
+    # create a som map
+    som_map = sits_som_map(samples_modis_ndvi)
+    # plot the som map
+    plot(som_map)
+    # evaluate the som map and create clusters
+    clusters_som = sits_som_evaluate_cluster(som_map)
+    # plot the cluster evaluation
+    plot(clusters_som)
+    # clean the samples
+    new_samples = sits_som_clean_samples(som_map)
diff --git a/pysits/docs/content/sits_som_evaluate_cluster.md b/pysits/docs/content/sits_som_evaluate_cluster.md
new file mode 100644
index 0000000..da98435
--- /dev/null
+++ b/pysits/docs/content/sits_som_evaluate_cluster.md
@@ -0,0 +1,26 @@
+Evaluate cluster
+
+`sits_som_evaluate_cluster()` produces a `SITSFrame` with the clusters
+found by the SOM map. For each cluster, it provides the percentage of
+classes inside it.
+
+Args:
+    som_map (SITSFrame): A SOM map produced by the `sits_som_map()`
+        function.
+
+Returns:
+    SITSFrame: The purity for each cluster.
+
+Examples:
+    from pysits import *
+
+    # create a som map
+    som_map = sits_som_map(samples_modis_ndvi)
+    # plot the som map
+    plot(som_map)
+    # evaluate the som map and create clusters
+    clusters_som = sits_som_evaluate_cluster(som_map)
+    # plot the cluster evaluation
+    plot(clusters_som)
+    # clean the samples
+    new_samples = sits_som_clean_samples(som_map)
diff --git a/pysits/docs/content/sits_som_map.md b/pysits/docs/content/sits_som_map.md
new file mode 100644
index 0000000..edf0f2f
--- /dev/null
+++ b/pysits/docs/content/sits_som_map.md
@@ -0,0 +1,76 @@
+Build a SOM for quality analysis of time series samples
+
+These function use self-organized maps to perform quality analysis in satellite
+image time series.
+
+Args:
+    data (SITSTimeSeriesModel): samples to be clustered.
+    grid_xdim (int): X dimension of the SOM grid (default = 25).
+    grid_ydim (int): Y dimension of the SOM grid.
+    alpha (float): Starting learning rate (decreases according to number
+        of iterations).
+    rlen (int): Number of iterations to produce the SOM.
+    distance (str): The type of similarity measure (distance). The
+        following similarity measurements are supported: `"euclidean"`,
+        `"dtw"`, and `"cosine"`. The default similarity measure is
+        `"dtw"`. For single-timestep samples (e.g. embeddings), all bands
+        are treated as one feature vector. `"dtw"` is not applicable and
+        falls back to `"euclidean"`, while `"cosine"` compares the angle
+        between the full embedding vectors.
+    som_radius (float): Radius of SOM neighborhood.
+    mode (str): Type of learning algorithm. The following learning
+        algorithm are available: `"online"`, `"batch"`, and `"pbatch"`.
+        The default learning algorithm is `"online"`.
+
+Returns:
+    SITStructureData: a structure with three members: (1) the samples,
+    with one additional column indicating to which neuron each sample has
+    been mapped; (2) the Kohonen map, used for plotting and cluster
+    quality measures; (3) the labelled neurons, where each class of each
+    neuron is associated to two values: (a) the prior probability that
+    this class belongs to a cluster based on the frequency of samples of
+    this class allocated to the neuron; (b) the posterior probability
+    that this class belongs to a cluster, using data for the neighbours
+    on the SOM map.
+
+Notes:
+    `sits_som_map` creates a SOM map, where high-dimensional data is mapped
+    into a two dimensional map, keeping the topological relations between data
+    patterns. Each sample is assigned to a neuron, and neurons are placed in
+    the grid based on similarity.
+    `sits_som_evaluate_cluster` analyses the neurons of the SOM map, and builds
+    clusters based on them. Each cluster is a neuron or a set of neuron
+    categorized with same label. It produces the percentage of mixture of
+    classes in each cluster.
+    `sits_som_clean_samples` evaluates sample quality based on the results of
+    the SOM map. The algorithm identifies noisy samples, using
+    `prior_threshold` for the prior probability and `posterior_threshold` for
+    the posterior probability. Each sample receives an evaluation tag,
+    according to the following rule: (a) If the prior probability is <
+    `prior_threshold`, the sample is tagged as "remove"; (b) If the prior
+    probability is >= `prior_threshold` and the posterior probability is
+    >=`posterior_threshold`, the sample is tagged as "clean"; (c) If the prior
+    probability is >= `posterior_threshold` and the posterior probability is <
+    `posterior_threshold`, the sample is tagged as "analyze" for further
+    inspection.
+    The user can define which tagged samples will be returned using the "keep"
+    parameter, with the following options: "clean", "analyze", "remove".
+    To learn more about the learning algorithms, check the `kohonen::supersom`
+    function.
+    The `sits` package implements the `"dtw"` (Dynamic Time Warping) similarity
+    measure. The `"euclidean"` similarity measurement come from the
+    `kohonen::supersom (dist.fcts)` function.
+
+Examples:
+    from pysits import *
+
+    # create a som map
+    som_map = sits_som_map(samples_modis_ndvi)
+    # plot the som map
+    plot(som_map)
+    # evaluate the som map and create clusters
+    clusters_som = sits_som_evaluate_cluster(som_map)
+    # plot the cluster evaluation
+    plot(clusters_som)
+    # clean the samples
+    new_samples = sits_som_clean_samples(som_map)
diff --git a/pysits/docs/content/sits_ssl_lejepa.md b/pysits/docs/content/sits_ssl_lejepa.md
new file mode 100644
index 0000000..116951a
--- /dev/null
+++ b/pysits/docs/content/sits_ssl_lejepa.md
@@ -0,0 +1,76 @@
+Self-supervised LeJEPA pre-training with resampling augmentation
+
+Self-supervised pre-training using the LeJEPA (Lean Joint-Embedding Predictive
+Architecture) loss and a torch encoder. Two views of each sample are created
+on-the-fly using the resampling augmentation of Saget et al. (2025). No labels
+are required.
+The LeJEPA loss combines two terms:
+1. Invariance: MSE between each view's projected embedding and the mean across
+   views \u2014 pulls views together.
+2. SIGReg (Sketched Isotropic Gaussian Regularization): constrains embeddings
+   to an isotropic Gaussian distribution by comparing 1-D projections to the
+   Gaussian characteristic function \u2014 prevents representational collapse
+   without teacher-student networks or stop-gradients.
+These two terms are balanced by a single hyperparameter `lambda`.
+The function can be used in two ways:
+- If `samples` is provided, it trains immediately and returns an encoder-ready
+  model object (see Value).
+- If `samples = None`, it returns a training function with signature
+  `function(samples)` that can be passed to `sits_pre_train` or called later.
+
+The augmentation strategy uses the resampling method of Saget et al. (2025):
+the time series is upsampled, two disjoint subsequences are drawn with a
+temporal coverage constraint, and each is resampled back to the original
+length.
+Unlike VICReg (which uses three separate regularization terms), LeJEPA replaces
+all collapse-prevention heuristics with the single SIGReg objective, yielding a
+simpler and more theoretically grounded method with a single trade-off
+hyperparameter.
+
+Args:
+    samples (SITSTimeSeriesModel): Samples object. If `None` (default),
+        returns a training function. If provided, triggers immediate
+        training.
+    embedding_dim (int): Dimensionality of the encoder embedding. Default:
+        64.
+    proj_dim (int): Dimensionality of the projector head used only during
+        pre-training. Default: 128.
+    lambda (float): Value in (0, 1). Trade-off between invariance
+        (`1 - lambda`) and SIGReg (`lambda`). Default: 0.02.
+    num_knots (int): Number of quadrature knots for the SIGReg
+        characteristic function test. Default: 17.
+    num_slices (int): Number of random projection directions for SIGReg.
+        Default: 256.
+    encoder_model (SITSMachineLearningMethod): Deep learning method that
+        takes time series as input and produces latent representations that
+        are used to compute the loss function (suggested options:
+        `sits_tempcnn()`, `sits_lighttae()`, `sits_resnet()`). Default:
+        `sits_tempcnn()`.
+    epochs (int): Maximum number of training epochs.
+    batch_size (int): Batch size for training. Default: 512.
+    validation_split (float): Value in (0, 1). Fraction of samples held out
+        for validation loss monitoring.
+    optimizer: A `torch` optimizer constructor (default:
+        `torch::optim_adamw`).
+    opt_hparams (dict): Optimizer hyperparameters.
+    lr_decay_epochs (int): Step size (in epochs) for LR decay.
+    lr_decay_rate (float): Multiplicative LR decay factor.
+    patience (int): Early-stopping patience.
+    min_delta (float): Minimum improvement for early stopping.
+    verbose (bool): Print training progress?
+    seed (int): Random seed for reproducibility.
+
+Returns:
+    R: If `samples = None`, a training function. If `samples` is provided, a
+        pretrained encoder.
+
+Examples:
+    from pysits import *
+
+    model = sits_pre_train(
+        samples_modis_ndvi,
+        sits_ssl_lejepa(
+            embedding_dim=32,
+            epochs=20
+        )
+    )
diff --git a/pysits/docs/content/sits_ssl_mae.md b/pysits/docs/content/sits_ssl_mae.md
new file mode 100644
index 0000000..7d772e0
--- /dev/null
+++ b/pysits/docs/content/sits_ssl_mae.md
@@ -0,0 +1,91 @@
+Pre-train a Masked Autoencoder for self-supervised learning of time series
+
+`sits_ssl_mae()` Implements a masked autoencoder (MAE) algorithm. It performs
+self-supervised learning by masking a subset of timesteps in each sample time
+series, training an encoder-decoder model to reconstruct the original signal,
+and returning the pretrained encoder as a `torch` module.
+The function can be used in two ways:
+- If `samples` is provided, it trains immediately and returns an encoder-ready
+  model object (see Value).
+- If `samples = None`, it returns a training function with signature
+  `function(samples)` that can be passed to `sits_pre_train` or called later.
+
+During training, the procedure:
+1. Masks each sample time series according to `masking_method`, `mask_ratio`,
+   `mask_value`, and `masked_bands`.
+2. Normalizes inputs and targets using quantile-based statistics derived from
+   the original samples.
+3. Splits data into training and validation partitions according to
+   `validation_split`.
+4. Trains an encoder-decoder model with a masked MSE objective, where loss is
+   computed only over masked positions.
+5. Applies early stopping and step learning-rate scheduling.
+6. Discards the decoder after training, returning the pretrained encoder for
+   use in downstream tasks.
+The decoder used during pretraining is a multilayer perceptron (MLP)
+specifically designed for MAE reconstruction. This decoder maps latent
+embeddings back to full time-series representations before being discarded
+after pretraining.
+Note: unlike the original vision MAE (He et al., 2022), where the encoder
+processes only visible (unmasked) patches, this implementation feeds the full
+time series to the encoder with masked positions replaced by `mask_value`. This
+simplification is standard for temporal data and functions as a denoising
+autoencoder variant of the MAE objective.
+When GPU execution is enabled in the environment, training may run on GPU via
+`luz` accelerators. Otherwise, it runs on CPU.
+
+Args:
+    samples (SITSTimeSeriesModel): Sample time series. If `None` (default),
+        returns a training function. If provided, triggers immediate
+        training. Base data samples (e.g., `sits_base`) are not supported.
+    embedding_dim (int): Dimensionality of the latent embedding produced by
+        the encoder (Default: 32).
+    encoder_model (SITSMachineLearningMethod): Deep learning method that takes
+        time series as input and produces latent representations that are used
+        to compute the loss function (suggested options: `sits_tempcnn()`,
+        `sits_lighttae()`, `sits_resnet()`). Default: `sits_tempcnn()`.
+    decoder_width (int): Width of the decoder MLP hidden layer.
+    dropout_rate (float): Dropout rates (0,1) for the linear module of the
+        decoder.
+    masking_method (str): Mask selection strategy. Options are `"random"` or
+        `"contiguous"`.
+    mask_ratio (float): Fraction of timesteps to mask, in (0, 1).
+    mask_value (float): Fill value used for masked timesteps.
+    masked_bands (list[str]): Which bands to mask. If `None`, all bands are
+        eligible for masking.
+    epochs (int): Maximum number of training epochs.
+    batch_size (int): Batch size used for training and validation.
+    validation_split (float): Fraction of samples held out for validation loss
+        monitoring, in (0, 1).
+    optimizer: A `torch` optimizer constructor, such as `torch::optim_adamw`.
+    opt_hparams (dict): Optimizer hyperparameters passed to `optimizer`.
+        Common entries include `lr`, `eps`, and `weight_decay`. Only
+        parameters supported by the chosen optimizer are accepted.
+    lr_decay_epochs (int): Step size (in epochs) for learning-rate decay when
+        using the step scheduler.
+    lr_decay_rate (float): Multiplicative decay factor applied by the
+        learning-rate scheduler.
+    patience (int): Number of epochs without improvement in validation loss
+        before early stopping.
+    min_delta (float): Minimum decrease in validation loss required to reset
+        the early-stopping patience counter.
+    verbose (bool): If `True`, prints training progress and per-epoch losses.
+    seed (int): Random seed used to initialize Torch randomness.
+
+Returns:
+    R: If `samples = None`, a training function with signature
+    `function(samples)` that trains an MAE and returns a pretrained encoder (a
+    `torch` module). If `samples` is provided, the result of applying the
+    training function to `samples` (i.e., a pretrained encoder-ready model
+    object used by the pretraining pipeline).
+
+Examples:
+    from pysits import *
+
+    model = sits_pre_train(
+        samples_modis_ndvi,
+        sits_ssl_mae(
+            embedding_dim=32,
+            epochs=20
+        )
+    )
diff --git a/pysits/docs/content/sits_ssl_vicreg.md b/pysits/docs/content/sits_ssl_vicreg.md
new file mode 100644
index 0000000..b1b785d
--- /dev/null
+++ b/pysits/docs/content/sits_ssl_vicreg.md
@@ -0,0 +1,87 @@
+Self-supervised VICReg pre-training with time-warping augmentation
+
+Self-supervised pre-training using the VICReg (Variance-Invariance-Covariance
+Regularization) loss and a torch encoder. Two views of each sample are created
+on-the-fly using the resampling augmentation of Saget et al. (2025): the time
+series is upsampled, two disjoint subsequences are drawn with a temporal
+coverage constraint, and each is resampled back to the original length. No
+labels are required.
+Both views are passed through a shared encoder and projector. The VICReg loss
+combines three objectives:
+1. Invariance: MSE between the projected representations of the two views 
+
+ce loss combines three objectives:
+1. Invariance: MSE between the projected representations of the two views 
+   pulls paired embeddings together.
+2. Variance: a hinge loss that keeps the standard deviation of each embedding
+   feature above a threshold of 1 across the batch 
+   prevents informational collapse.
+3. Covariance: penalizes off-diagonal entries of the embedding covariance
+   matrix 
+   decorrelates features.
+The function can be used in two ways:
+- If `samples` is provided, it trains immediately and returns an encoder-ready
+  model object (see Returns).
+- If `samples` is `None`, it returns a training function that can be passed to
+  `sits_pre_train` or called later.
+
+The augmentation strategy creates two views of each sample using the resampling
+method of Saget et al. (2025). The original time series (length `T`) is
+upsampled to `2T` timesteps by linear interpolation. Two disjoint subsequences
+of `T/2` timesteps are drawn from the upsampled series, with a constraint that
+at least `floor((T/2)/4)` timesteps fall in each temporal quarter. Each
+subsequence is then resampled back to `T` positions by rescaling its timestamps
+and interpolating, producing two views that preserve overall temporal structure
+while differing in fine-grained detail. Because the subsampling is random,
+views differ at every epoch.
+
+Args:
+    samples (SITSTimeSeriesModel): Samples object. If `None` (default),
+        returns a training function. If provided, triggers immediate
+        training. Base data samples are not supported.
+    embedding_dim (int): Dimensionality of the encoder embedding (exported
+        features). Default: 64.
+    proj_dim (int): Dimensionality of the projector head used only during
+        pre-training. Default: 128.
+    sim_coeff (float): Weight of the invariance (MSE) term in the VICReg
+        loss. Default: 25.0.
+    std_coeff (float): Weight of the variance (hinge) term in the VICReg
+        loss. Default: 25.0.
+    cov_coeff (float): Weight of the covariance (off-diagonal) term in the
+        VICReg loss. Default: 1.0.
+    encoder_model (SITSMachineLearningMethod): Deep learning method that
+        takes time series as input and produces latent representations that
+        are used to compute the loss function (suggested options:
+        `sits_tempcnn()`, `sits_lighttae()`, `sits_resnet()`). Default:
+        `sits_tempcnn()`.
+    epochs (int): Maximum number of training epochs.
+    batch_size (int): Batch size for training. Default: 128.
+    validation_split (float): Fraction of samples held out for validation
+        loss monitoring, in (0, 1).
+    optimizer: A `torch` optimizer constructor (default:
+        `torch::optim_adamw`).
+    opt_hparams (dict): Optimizer hyperparameters. Common entries: `lr`,
+        `eps`, `weight_decay`.
+    lr_decay_epochs (int): Step size (in epochs) for LR decay.
+    lr_decay_rate (float): Multiplicative LR decay factor.
+    patience (int): Early-stopping patience (epochs without improvement).
+    min_delta (float): Minimum improvement required to reset the patience
+        counter.
+    verbose (bool): Print training progress?
+    seed (int): Random seed for reproducibility.
+
+Returns:
+    R: If `samples` is `None`, a training function that trains a VICReg model
+    and returns a pretrained encoder. If `samples` is provided, the result of
+    applying the training function to `samples` directly.
+
+Examples:
+    from pysits import *
+
+    model = sits_pre_train(
+        samples_modis_ndvi,
+        sits_ssl_vicreg(
+            embedding_dim=32,
+            epochs=20
+        )
+    )
diff --git a/pysits/docs/content/sits_stats.md b/pysits/docs/content/sits_stats.md
new file mode 100644
index 0000000..68fae17
--- /dev/null
+++ b/pysits/docs/content/sits_stats.md
@@ -0,0 +1,21 @@
+Obtain statistics for all sample bands
+
+Most machine learning algorithms require data to be normalized. This
+applies to the "SVM" method and to all deep learning ones. To
+normalize the predictors, it is necessary to extract the statistics
+of each band of the samples. This function computes the 2 of the
+distribution of each band of the samples. This values are used as
+minimum and maximum values in the normalization operation performed
+by the sits_pred_normalize() function.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series samples used as
+        training data.
+
+Returns:
+    SITStructureData: The training data statistics.
+
+Examples:
+    from pysits import *
+
+    stats = sits_stats(samples_modis_ndvi)
diff --git a/pysits/docs/content/sits_stratified_sampling.md b/pysits/docs/content/sits_stratified_sampling.md
new file mode 100644
index 0000000..1028f12
--- /dev/null
+++ b/pysits/docs/content/sits_stratified_sampling.md
@@ -0,0 +1,73 @@
+Retrieval of sample locations by strata for a classified cube
+
+This function can be used in two ways: (a) When the parameter
+"sampling_design" is available, it takes the cube with different labels
+and a column for the sampling design table with a number of samples per
+class and allocates a set of locations for each class.
+(b) When the parameter "sampling_design" is not provided, the method
+selects a set of locations based on a classified cube based on the
+parameter "samples_per_class". This parameter should either be a `dict`
+(keys are the labels of the cube) or a single value. In the latter case,
+the same value is used to retrieve the samples for all classes.
+
+Args:
+    cube (SITSCubeModel): Classified cube.
+    sampling_design (pandas.DataFrame): Result of `sits_sampling_design`.
+    alloc (str): Allocation method chosen.
+    samples_per_class (int | dict): Number of samples per class (in case
+        `sampling_design` is `None`). Either a single integer (applied to
+        all classes) or a `dict` keyed by class label. When keyed, keys
+        must be valid class labels of the cube but do not need to cover
+        all classes — only the listed classes will be sampled.
+    overhead (float): Additional percentage to account for border points.
+    multicores (int): Number of cores that will be used to sample the
+        images in parallel.
+    memsize (int): Memory available for sampling.
+    shp_file (str | pathlib.Path): Name of shapefile to be saved
+        (optional).
+    progress (bool): Show progress bar? Default is `True`.
+
+Returns:
+    SITSFrameSF: Point object with required samples and label.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # label the probability cube
+    label_cube = sits_label_classification(
+        probs_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # Option 1 - select samples based on sampling design
+    # estimated UA for classes
+    expected_ua = {
+        "Cerrado": 0.95, "Forest": 0.95,
+        "Pasture": 0.95, "Soy_Corn": 0.95
+    }
+    # design sampling
+    sampling_design = sits_sampling_design(label_cube, expected_ua)
+    # select samples using the sampling design
+    samples = sits_stratified_sampling(
+        label_cube,
+        sampling_design=sampling_design,
+        alloc="alloc_prop"
+    )
+    # Option 2 - Select samples based on a fixed number of samples per class
+    samples = sits_stratified_sampling(
+        label_cube,
+        samples_per_class=100
+    )
diff --git a/pysits/docs/content/sits_svm.md b/pysits/docs/content/sits_svm.md
new file mode 100644
index 0000000..d15839a
--- /dev/null
+++ b/pysits/docs/content/sits_svm.md
@@ -0,0 +1,52 @@
+Train support vector machine models
+
+This function receives a set of training samples with a set of
+attributes X for each observation Y. These attributes are the values of
+the time series for each band. The SVM algorithm is used for
+multiclass-classification. For this purpose, it uses the
+"one-against-one" approach, in which k(k-1)/2 binary classifiers are
+trained; the appropriate class is found by a voting scheme. This
+function is a front-end to the "svm" method in the "e1071" package.
+Please refer to the documentation in that package for more details.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    formula: Symbolic description of the model to be fit. (default:
+        sits_formula_linear).
+    scale (list[bool]): Indicates the variables to be scaled.
+    cachesize (float): Cache memory in MB (default = 1000).
+    kernel (str): Kernel used in training and predicting. options:
+        "linear", "polynomial", "radial", "sigmoid" (default: "radial").
+    degree (int): Exponential of polynomial type kernel (default: 3).
+    coef0 (float): Parameter needed for kernels of type polynomial and
+        sigmoid (default: 0).
+    cost (float): Cost of constraints violation (default: 10).
+    tolerance (float): Tolerance of termination criterion
+        (default: 0.001).
+    epsilon (float): Epsilon in the insensitive-loss function
+        (default: 0.1).
+    cross (int): Number of cross validation folds applied to assess the
+        quality of the model (default: 10).
+    **kwargs (dict): Other parameters to be passed to e1071::svm
+        function.
+
+Returns:
+    SITSMachineLearningMethod: Model fitted to input data (to be passed
+        to `sits_classify`).
+
+Notes:
+    Please refer to the sits documentation available in
+    https://e-sensing.github.io/sitsbook/ for detailed examples.
+
+Examples:
+    from pysits import *
+
+    # Example of training a model for time series classification
+    # Retrieve the samples for Mato Grosso
+    # train an SVM model
+    ml_model = sits_train(samples_modis_ndvi, ml_method=sits_svm)
+    # classify the point
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    # classify the point
+    point_class = sits_classify(data=point_ndvi, ml_model=ml_model)
+    plot(point_class)
diff --git a/pysits/docs/content/sits_tae.md b/pysits/docs/content/sits_tae.md
new file mode 100644
index 0000000..b384fe6
--- /dev/null
+++ b/pysits/docs/content/sits_tae.md
@@ -0,0 +1,80 @@
+Train a model using Temporal Self-Attention Encoder
+
+Implementation of Temporal Attention Encoder (TAE) for satellite image time
+series classification.
+TAE is a simplified version of the well-known self-attention architecture used
+in large language models. Its modified self-attention scheme that uses the
+input embeddings as values. TAE defines a single master query for each
+sequence, computed from the temporal average of the queries. This master query
+is compared to the sequence of keys to produce a single attention mask used to
+weight the temporal mean of values into a single feature vector.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    samples_validation (SITSTimeSeriesModel): Time series with the
+        validation samples. If the `samples_validation` parameter is
+        provided, the `validation_split` parameter is ignored.
+    epochs (int): Number of iterations to train the model.
+    batch_size (int): Number of samples per gradient update.
+    validation_split (float): Number between 0 and 1. Fraction of training
+        data to be used as validation data.
+    optimizer: Optimizer function to be used.
+    opt_hparams (dict): Hyperparameters for optimizer: lr : Learning rate of
+        the optimizer eps: Term added to the denominator to improve
+        numerical stability. weight_decay: L2 regularization
+    lr_decay_epochs (int): Number of epochs to reduce learning rate.
+    lr_decay_rate (float): Decay factor for reducing learning rate.
+    patience (int): Number of epochs without improvements until training
+        stops.
+    min_delta (float): Minimum improvement to reset the patience counter.
+    seed (int): Seed for random values.
+    verbose (bool): Verbosity mode. Default is False.
+
+Returns:
+    SITSMachineLearningMethod: A fitted model to be used for classification.
+
+Notes:
+    `sits` provides a set of default values for all classification models.
+    These settings have been chosen based on testing by the authors.
+    Nevertheless, users can control all parameters for each model. Novice users
+    can rely on the default values, while experienced ones can fine-tune deep
+    learning models using `sits_tuning`.
+    This function is based on the paper by Vivien Garnot referenced below and
+    code available on github at https://github.com/VSainteuf/pytorch-psetae.
+    We also used the code made available by Maja Schneider in her work with
+    Marco K
+    https://github.com/maja601/RC2020-psetae.
+    If you use this method, please cite Garnot's and Schneider's work.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a TAE model
+    torch_model = sits_train(samples_modis_ndvi, sits_tae())
+    # plot the model
+    plot(torch_model)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=torch_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_tempcnn.md b/pysits/docs/content/sits_tempcnn.md
new file mode 100644
index 0000000..27e92df
--- /dev/null
+++ b/pysits/docs/content/sits_tempcnn.md
@@ -0,0 +1,83 @@
+Train temporal convolutional neural network models
+
+Use a TempCNN algorithm to classify data, which has two stages: a 1D CNN and a
+multi-layer perceptron. Users can define the depth of the 1D network, as well
+as the number of perceptron layers.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    samples_validation (SITSTimeSeriesModel): Time series with the validation
+        samples. If provided, the `validation_split` parameter is ignored.
+    cnn_layers (list[int]): Number of 1D convolutional filters per layer.
+    cnn_kernels (list[int]): Size of the 1D convolutional kernels.
+    cnn_dropout_rates (list[float]): Dropout rates for 1D convolutional
+        filters.
+    dense_layer_nodes (int): Number of nodes in the dense layer.
+    dense_layer_dropout_rate (float): Dropout rate (0,1) for the dense layer.
+    epochs (int): Number of iterations to train the model.
+    batch_size (int): Number of samples per gradient update.
+    validation_split (float): Fraction of training data to be used for
+        validation.
+    optimizer: Optimizer function to be used.
+    opt_hparams (dict): Hyperparameters for optimizer: lr : Learning rate of
+        the optimizer eps: Term added to the denominator to improve numerical
+        stability. weight_decay: L2 regularization
+    lr_decay_epochs (int): Number of epochs to reduce learning rate.
+    lr_decay_rate (float): Decay factor for reducing learning rate.
+    patience (int): Number of epochs without improvements until training stops.
+    min_delta (float): Minimum improvement in loss function to reset the
+        patience counter.
+    seed (int): Seed for random values.
+    verbose (bool): Verbosity mode. Default is `False`.
+
+Returns:
+    R: A fitted model to be used for classification.
+
+Notes:
+    `sits` provides a set of default values for all classification models.
+    These settings have been chosen based on testing by the authors.
+    Nevertheless, users can control all parameters for each model. Novice users
+    can rely on the default values, while experienced ones can fine-tune deep
+    learning models using `sits_tuning`.
+    This function is based on the paper by Charlotte Pelletier referenced
+    below. If you use this method, please cite the original tempCNN paper.
+    The torch version is based on the code made available by the BreizhCrops
+    team: Marc Russwurm, Charlotte Pelletier, Marco Korner, Maximilian Zollner.
+    The original python code is available at the website
+    https://github.com/dl4sits/BreizhCrops. This code is licensed as GPL-3.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a TempCNN model
+    torch_model = sits_train(
+        samples_modis_ndvi,
+        sits_tempcnn(epochs=20, verbose=True)
+    )
+    # plot the model
+    plot(torch_model)
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=torch_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    bayes_cube = sits_smooth(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the smoothed cube
+    plot(bayes_cube)
+    # label the probability cube
+    label_cube = sits_label_classification(
+        bayes_cube,
+        output_dir=tempfile.gettempdir()
+    )
+    # plot the labelled cube
+    plot(label_cube)
diff --git a/pysits/docs/content/sits_texture.md b/pysits/docs/content/sits_texture.md
new file mode 100644
index 0000000..02a2b99
--- /dev/null
+++ b/pysits/docs/content/sits_texture.md
@@ -0,0 +1,52 @@
+Apply a set of texture measures on a data cube.
+
+A set of texture measures based on the Grey Level Co-occurrence Matrix
+(GLCM) described by Haralick. Our implementation follows the guidelines
+and equations described by Hall-Beyer (both are referenced below).
+
+The spatial relation between the central pixel and its neighbor is
+expressed in radians values, where: #'
+- `0`: corresponds to the neighbor on right-side
+- `pi/4`: corresponds to the neighbor on the top-right diagonals
+- `pi/2`: corresponds to the neighbor on above
+- `3*pi/4`: corresponds to the neighbor on the top-left diagonals
+Our implementation relies on a symmetric co-occurrence matrix, which
+considers the opposite directions of an angle. For example, the neighbor
+pixels based on `0` angle rely on the left and right direction; the
+neighbor pixels of `pi/2` are above and below the central pixel, and so
+on. If more than one angle is provided, we compute their average.
+
+Args:
+    cube (SITSCubeModel): Valid data cube.
+    window_size (int): An odd number representing the size of the sliding
+        window.
+    angles (list[float]): The direction angles in radians related to the
+        central pixel and its neighbor (See details). Default is 0.
+    memsize (int): Memory available for classification (in GB).
+    multicores (int): Number of cores to be used for classification.
+    output_dir (str | pathlib.Path): Directory where files will be saved.
+    progress (bool): Show progress bar?
+    **kwargs (dict): GLCM function (see details).
+
+Returns:
+    SITSCubeModel: A data cube with new bands, produced according to the
+        requested measure.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # Compute the NDVI variance
+    cube_texture = sits_texture(
+        cube=cube,
+        NDVIVAR="glcm_variance(NDVI)",
+        window_size=5,
+        output_dir=tempfile.mkdtemp()
+    )
diff --git a/pysits/docs/content/sits_tiles_to_roi.md b/pysits/docs/content/sits_tiles_to_roi.md
new file mode 100644
index 0000000..acb4398
--- /dev/null
+++ b/pysits/docs/content/sits_tiles_to_roi.md
@@ -0,0 +1,25 @@
+Convert tile information to ROI in WGS84
+
+Takes a list of tiles from a given grid system and produces a ROI
+(region of interest) in WGS84 covering them.
+
+Args:
+    tiles (list[str]): Names of tiles from the selected `grid_system`.
+    grid_system (str): Grid system that the `tiles` belong to.
+        Currently supported grid systems are the MGRS grid
+        (`"MGRS"`, default) and those used by the Brazil Data Cube
+        (`"BDC_LG_V2"`, `"BDC_MD_V2"` and `"BDC_SM_V2"`).
+
+Returns:
+    SITSNamedVector: Valid ROI to use in other SITS functions.
+
+Examples:
+    from pysits import *
+
+    # Convert MGRS tiles to a ROI
+    roi = sits_tiles_to_roi(["22KGA", "22KGV"])
+
+    # Convert Brazil Data Cube (large grid) tiles to a ROI
+    roi = sits_tiles_to_roi(
+        ["003004", "003005"], grid_system="BDC_LG_V2"
+    )
diff --git a/pysits/docs/content/sits_timeline.md b/pysits/docs/content/sits_timeline.md
new file mode 100644
index 0000000..3f598e2
--- /dev/null
+++ b/pysits/docs/content/sits_timeline.md
@@ -0,0 +1,15 @@
+Get timeline of a cube or a set of time series
+
+This function returns the timeline for a given data set, either a set of
+time series, a data cube, or a trained model.
+
+Args:
+    data (SITSTimeSeriesModel | SITSCubeModel): time series or data cube.
+
+Returns:
+    list: timeline of samples or data cube.
+
+Examples:
+    from pysits import *
+
+    sits_timeline(samples_modis_ndvi)
diff --git a/pysits/docs/content/sits_to_csv.md b/pysits/docs/content/sits_to_csv.md
new file mode 100644
index 0000000..9f00b62
--- /dev/null
+++ b/pysits/docs/content/sits_to_csv.md
@@ -0,0 +1,23 @@
+Export sits time series metadata to the CSV format
+
+Converts metadata from a `SITSTimeSeriesModel` to a CSV file. The CSV
+file will not contain the actual time series. Its columns will be the
+same as those of a CSV file used to retrieve data from ground
+information ("latitude", "longitude", "start_date", "end_date",
+"cube", "label"). If the file is `None`, returns a `pandas.DataFrame`.
+
+Args:
+    data (SITSTimeSeriesModel): Time series.
+    file (str | pathlib.Path): Full path of the exported CSV file (valid
+        file name with extension ".csv").
+
+Returns:
+    SITSFrame: Data with CSV columns (optional).
+
+Examples:
+    from pysits import *
+    import tempfile
+    import os
+
+    csv_file = os.path.join(tempfile.gettempdir(), "cerrado_2classes.csv")
+    sits_to_csv(cerrado_2classes, file=csv_file)
diff --git a/pysits/docs/content/sits_to_xlsx.md b/pysits/docs/content/sits_to_xlsx.md
new file mode 100644
index 0000000..7031bac
--- /dev/null
+++ b/pysits/docs/content/sits_to_xlsx.md
@@ -0,0 +1,42 @@
+Save accuracy assessments as Excel files
+
+Saves confusion matrices as Excel spreadsheets. This function takes a
+list of accuracy assessments generated by `sits_accuracy` and saves
+them in an Excel spreadsheet.
+
+Args:
+    acc (SITSConfusionMatrix | list[SITSConfusionMatrix]): Accuracy
+        statistics, either an output of `sits_accuracy` or a list of
+        those.
+    file (str | pathlib.Path): File where the XLSX data is to be
+        saved.
+
+Returns:
+    None: Called for side effects.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # A dataset containing a tibble with time series samples
+    # for the Mato Grosso state in Brasil
+    # create a list to store the results
+    results = []
+
+    # accuracy assessment lightTAE
+    acc_ltae = sits_kfold_validate(samples_modis_ndvi,
+        folds=5,
+        multicores=1,
+        ml_method=sits_lighttae()
+    )
+    # use a name
+    acc_ltae["name"] = "LightTAE"
+
+    # put the result in a list
+    results.append(acc_ltae)
+
+    # save to xlsx file
+    sits_to_xlsx(
+        results,
+        file=tempfile.NamedTemporaryFile(prefix="accuracy_mato_grosso_dl_", suffix=".xlsx").name
+    )
diff --git a/pysits/docs/content/sits_train.md b/pysits/docs/content/sits_train.md
new file mode 100644
index 0000000..9c45e98
--- /dev/null
+++ b/pysits/docs/content/sits_train.md
@@ -0,0 +1,61 @@
+Train classification models
+
+Given a set of time series, returns trained models. Currently,
+sits supports the following models:
+- support vector machines: `sits_svm`;
+- random forests: `sits_rfor`;
+- extreme gradient boosting: `sits_xgboost`;
+- light gradient boosting: `sits_lightgbm`;
+- multi-layer perceptrons: `sits_mlp`;
+- temporal CNN: `sits_tempcnn`;
+- residual network encoders: `sits_resnet`;
+- LSTM with convolutional networks: `sits_lstm_fcn`;
+- temporal self-attention encoders: `sits_lighttae` and `sits_tae`.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    ml_method (SITSMachineLearningMethod): Machine learning method.
+
+Returns:
+    SITSMachineLearningMethod: Model fitted to input data to be passed to
+        `sits_classify`.
+
+Notes:
+    The main `sits` classification workflow has the following steps:
+    1. `sits_cube`: selects a ARD image collection from a cloud provider.
+    2. `sits_cube_copy`: copies an ARD image collection from a cloud provider
+       to a local directory for faster processing.
+    3. `sits_regularize`: create a regular data cube from an ARD image
+       collection.
+    4. `sits_apply`: create new indices by combining bands of a regular data
+       cube (optional).
+    5. `sits_get_data`: extract time series from a regular data cube based on
+       user-provided labelled samples.
+    6. `sits_train`: train a machine learning model based on image time series.
+    7. `sits_classify`: classify a data cube using a machine learning model and
+       obtain a probability cube.
+    8. `sits_smooth`: post-process a probability cube using a spatial smoother
+       to remove outliers and increase spatial consistency.
+    9. `sits_label_classification`: produce a classified map by selecting the
+       label with the highest probability from a smoothed cube.
+    `sits_train` provides a standard interface to machine learning models. It
+    takes two mandatory parameters: the training data (`samples`) and the ML
+    algorithm (`ml_method`). The output is a model that can be used to classify
+    individual time series or data cubes with `sits_classify`.
+    `sits` provides a set of default values for all classification models.
+    These settings have been chosen based on testing by the authors.
+    Nevertheless, users can control all parameters for each model. Novice users
+    can rely on the default values, while experienced ones can fine-tune deep
+    learning models using `sits_tuning`.
+
+Examples:
+    from pysits import *
+
+    # Retrieve the set of samples for Mato Grosso
+    # fit a training model (rfor model)
+    ml_model = sits_train(samples_modis_ndvi, sits_rfor(num_trees=50))
+    # get a point and classify the point with the ml_model
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    class_ = sits_classify(
+        data=point_ndvi, ml_model=ml_model
+    )
diff --git a/pysits/docs/content/sits_tuning.md b/pysits/docs/content/sits_tuning.md
new file mode 100644
index 0000000..25eabe8
--- /dev/null
+++ b/pysits/docs/content/sits_tuning.md
@@ -0,0 +1,61 @@
+Tuning machine learning models hyper-parameters
+
+This function performs a random search on values of selected hyperparameters,
+and produces a data frame with the accuracy and kappa values produced by a
+validation procedure. The result allows users to select appropriate
+hyperparameters for deep learning models.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series set to be validated.
+    samples_validation (SITSTimeSeriesModel): Time series set used for
+        validation.
+    validation_split (float): Percent of original time series set to be
+        used for validation (if `samples_validation` is `None`).
+    ml_method (SITSMachineLearningMethod): Machine learning method.
+    params (dict): Hyper parameters to be passed to `ml_method`. User can
+        use `uniform`, `choice`, `randint`, `normal`, `lognormal`,
+        `loguniform`, and `beta` distribution functions to randomize
+        parameters.
+    trials (int): Number of random trials to perform the search.
+    multicores (int): Number of cores to process in parallel.
+    gpu_memory (int): Memory available in GPU in GB (default = 4).
+    batch_size (int): Batch size for GPU classification.
+    progress (bool): Show progress bar?
+
+Returns:
+    SITSTuningResults: All parameters used to train on each trial ordered
+        by accuracy.
+
+Notes:
+    Machine learning algorithms have hyperparameters that control the
+    algorithm's behaviour. This function allows users to test different
+    combinations of hyperparameters for a given sample set, thus selecting a
+    set of values which fits the training data. The `sits_tuning` function can
+    be used with both traditional machine learning methods (e.g., random
+    forests) as well as deep learning ones.
+    Instead of performing an exhaustive test of all parameter combinations,
+    `sits_tuning` selects them randomly. Validation is done using an
+    independent set of samples or by a validation split. The function returns
+    the best hyper-parameters in a `dict`. Hyper-parameters passed to `params`
+    parameter should be passed by calling `sits_tuning_hparams`.
+    Deep learning models use stochastic gradient descent (SGD) techniques to
+    find optimal solutions. To perform SGD, models use optimization algorithms
+    which have hyperparameters that have to be adjusted to achieve best
+    performance for each application.
+    When using a GPU for deep learning, `gpu_memory` indicates the memory of
+    the graphics card which is available for processing. The parameter
+    `batch_size` defines the size of the matrix (measured in number of rows)
+    which is sent to the GPU for classification. Users can test different
+    values of `batch_size` to find out which one best fits their GPU
+    architecture.
+    It is not possible to have an exact idea of the size of Deep Learning
+    models in GPU memory, as the complexity of the model and factors such as
+    CUDA Context increase the size of the model in memory. Therefore, we
+    recommend that you leave at least 1GB free on the video card to store the
+    Deep Learning model that will be used.
+    For users of Apple M3 chips or similar with a Neural Engine, be aware that
+    these chips share memory between the GPU and the CPU. Tests indicate that
+    the `memsize` should be set to half to the total memory and the
+    `batch_size` parameter should be a small number (we suggest the value of
+    64). Be aware that increasing these parameters may lead to memory
+    conflicts.
diff --git a/pysits/docs/content/sits_tuning_hparams.md b/pysits/docs/content/sits_tuning_hparams.md
new file mode 100644
index 0000000..e184ae0
--- /dev/null
+++ b/pysits/docs/content/sits_tuning_hparams.md
@@ -0,0 +1,29 @@
+Tuning machine learning models hyper-parameters
+
+This function allow user building the hyper-parameters space used by
+`sits_tuning()` function search randomly the best parameter combination.
+Users should pass the possible values for hyper-parameters as constants or by
+calling the following random functions:
+- `uniform(min = 0, max = 1, n = 1)`: returns random numbers from a uniform
+  distribution with parameters min and max.
+- `choice(..., replace = TRUE, n = 1)`: returns random objects passed to `...`
+  with replacement or not (parameter `replace`).
+- `randint(min, max, n = 1)`: returns random integers from a uniform
+  distribution with parameters min and max.
+- `normal(mean = 0, sd = 1, n = 1)`: returns random numbers from a normal
+  distribution with parameters min and max.
+- `lognormal(meanlog = 0, sdlog = 1, n = 1)`: returns random numbers from a
+  lognormal distribution with parameters min and max.
+- `loguniform(minlog = 0, maxlog = 1, n = 1)`: returns random numbers from a
+  loguniform distribution with parameters min and max.
+- `beta(shape1, shape2, n = 1)`: returns random numbers from a beta
+  distribution with parameters min and max.
+These functions accepts `n` parameter to indicate how many values should be
+returned.
+
+Args:
+    **kwargs (dict): Used to prepare hyper-parameter space.
+
+Returns:
+    TuningFunctionCall: the hyper-parameter space to be passed to
+        `sits_tuning()`'s `params` parameter.
diff --git a/pysits/docs/content/sits_uncertainty.md b/pysits/docs/content/sits_uncertainty.md
new file mode 100644
index 0000000..1cf78e3
--- /dev/null
+++ b/pysits/docs/content/sits_uncertainty.md
@@ -0,0 +1,58 @@
+Estimate classification uncertainty based on probs cube
+
+Calculate the uncertainty cube based on the probabilities produced by the
+classifier. Takes a `probability cube` as input and produces a `uncertainty
+cube`.
+
+Args:
+    cube (SITSCubeModel): Probability data cube.
+    type (str): Method to measure uncertainty. See details.
+    multicores (int): Number of cores to run the function.
+    memsize (int): Maximum overall memory (in GB) to run the function.
+    output_dir (str | pathlib.Path): Output directory for image files.
+    version (str): Version of resulting image (in the case of multiple
+        tests).
+    progress (bool): Check progress bar?
+    agg_method (str): Aggregation method for probabilities. One of "mean"
+        (default) or "median".
+    **kwargs (dict): Other parameters for specific functions.
+
+Returns:
+    SITSCubeModel: An uncertainty data cube.
+
+Notes:
+    The output of `sits_classify` and `sits_smooth` is a `probability cube`
+    containing the class probability for all pixels, which are generated by the
+    machine learning model. The `sits_uncertainty` function takes a
+    `probability cube` and produces a `uncertainty code` which contains a
+    measure of uncertainty for each pixel, based on the class probabilities.
+    The uncertainty measure is relevant in the context of active leaning, and
+    helps to increase the quantity and quality of training samples by providing
+    information about the confidence of the model.
+    The supported types of uncertainty are:
+    1. `entropy`: the difference between all predictions expressed a Shannon
+       measure of entropy.
+    2. `least`: the difference between 1.0 and most confident prediction.
+    3. `margin`: the difference between the two most confident predictions.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # calculate uncertainty
+    uncert_cube = sits_uncertainty(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the resulting uncertainty cube
+    plot(uncert_cube)
diff --git a/pysits/docs/content/sits_uncertainty_sampling.md b/pysits/docs/content/sits_uncertainty_sampling.md
new file mode 100644
index 0000000..b0e80f5
--- /dev/null
+++ b/pysits/docs/content/sits_uncertainty_sampling.md
@@ -0,0 +1,62 @@
+Suggest samples for enhancing classification accuracy
+
+Suggest samples for regions of high uncertainty as predicted by the model. The
+function selects data points that have confused an algorithm. These points
+don't have labels and need be manually labelled by experts and then used to
+increase the classification's training set.
+This function is best used in the following context:
+1. Select an initial set of samples.
+2. Train a machine learning model.
+3. Build a data cube and classify it using the model.
+4. Run a Bayesian smoothing in the resulting probability cube.
+5. Create an uncertainty cube.
+6. Perform uncertainty sampling.
+The Bayesian smoothing procedure will reduce the classification outliers and
+thus increase the likelihood that the resulting pixels with high uncertainty
+have meaningful information.
+
+Args:
+    uncert_cube (SITSCubeModel): An uncertainty cube. See
+        `sits_uncertainty`.
+    n (int): Number of suggested points to be sampled per tile.
+    min_uncert (float): Minimum uncertainty value to select a sample.
+    max_uncert (float): Maximum uncertainty value to select a sample.
+        Default is Inf (no upper limit).
+    sampling_window (int): Window size for collecting points (in pixels).
+        The minimum window size is 10.
+    multicores (int): Number of workers for parallel processing (min = 1,
+        max = 2048).
+    memsize (int): Maximum overall memory (in GB) to run the function.
+    progress (bool): Whether to show progress bars.
+
+Returns:
+    SITSFrame: Longitude and latitude in WGS84 with locations which have
+        high uncertainty and meet the minimum distance criteria.
+
+Examples:
+    from pysits import *
+
+    # create a data cube
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # build a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, ml_method=sits_rfor())
+    # classify the cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempdir()
+    )
+    # create an uncertainty cube
+    uncert_cube = sits_uncertainty(probs_cube,
+        type="entropy",
+        output_dir=tempdir()
+    )
+    # obtain a new set of samples for active learning
+    # the samples are located in uncertain places
+    new_samples = sits_uncertainty_sampling(
+        uncert_cube,
+        n=10, min_uncert=0.4
+    )
diff --git a/pysits/docs/content/sits_validate.md b/pysits/docs/content/sits_validate.md
new file mode 100644
index 0000000..e80ba5f
--- /dev/null
+++ b/pysits/docs/content/sits_validate.md
@@ -0,0 +1,61 @@
+Validate time series samples
+
+One round of cross-validation involves partitioning a sample of data into
+complementary subsets, performing the analysis on one subset (called the
+training set), and validating the analysis on the other subset (called the
+validation set or testing set).
+The function takes two arguments: a set of time series with a machine learning
+model and another set with validation samples. If the validation sample set is
+not provided, The sample dataset is split into two parts, as defined by the
+parameter validation_split. The accuracy is determined by the result of the
+validation test set.
+This function returns the confusion matrix, and Kappa values.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series to be validated.
+    samples_validation (SITSTimeSeriesModel): Optional time series used for
+        validation.
+    validation_split (float): Percent of original time series set to be
+        used for validation if `samples_validation` is `None`.
+    ml_method (SITSMachineLearningMethod): Machine learning method.
+    gpu_memory (int): Memory available in GPU in GB (default = 4).
+    batch_size (int): Batch size for GPU classification.
+
+Returns:
+    SITSConfusionMatrix: Confusion matrix to be used for validation
+        assessment.
+
+Notes:
+    When using a GPU for deep learning, `gpu_memory` indicates the memory of
+    the graphics card which is available for processing. The parameter
+    `batch_size` defines the size of the matrix (measured in number of rows)
+    which is sent to the GPU for classification. Users can test different
+    values of `batch_size` to find out which one best fits their GPU
+    architecture.
+    It is not possible to have an exact idea of the size of Deep Learning
+    models in GPU memory, as the complexity of the model and factors such as
+    CUDA Context increase the size of the model in memory. Therefore, we
+    recommend that you leave at least 1GB free on the video card to store the
+    Deep Learning model that will be used.
+    For users of Apple M3 chips or similar with a Neural Engine, be aware that
+    these chips share memory between the GPU and the CPU. Tests indicate that
+    the `memsize` should be set to half to the total memory and the
+    `batch_size` parameter should be a small number (we suggest the value of
+    64). Be aware that increasing these parameters may lead to memory
+    conflicts.
+
+Examples:
+    from pysits import *
+
+    samples = sits_sample(cerrado_2classes, frac=0.5)
+    samples_validation = sits_sample(cerrado_2classes, frac=0.5)
+    conf_matrix_1 = sits_validate(
+        samples=samples,
+        samples_validation=samples_validation,
+        ml_method=sits_rfor()
+    )
+    conf_matrix_2 = sits_validate(
+        samples=cerrado_2classes,
+        validation_split=0.2,
+        ml_method=sits_rfor()
+    )
diff --git a/pysits/docs/content/sits_variance.md b/pysits/docs/content/sits_variance.md
new file mode 100644
index 0000000..abfd7c0
--- /dev/null
+++ b/pysits/docs/content/sits_variance.md
@@ -0,0 +1,50 @@
+Calculate the variance of a probability cube
+
+Takes a probability data cube (either a raster or a segmented
+vector cube) and estimates the variance of the logit of the
+probability. For a standard raster cube, this is a local sliding-window
+variance. For a vector/segmented cube, it calculates the variance of all
+pixels inside each segment. This supports the choice of parameters for
+Bayesian smoothing.
+
+Args:
+    cube (SITSCubeModel): Probability data cube.
+    window_size (int): Size of the neighborhood (odd integer). Not used
+        for a segmented vector cube.
+    neigh_fraction (float): Fraction of neighbors with highest
+        probability for Bayesian inference (from 0.0 to 1.0).
+    memsize (int): Maximum overall memory (in GB) to run the smoothing
+        (min = 1, max = 16384).
+    multicores (int): Number of cores to run the smoothing function
+        (min = 1, max = 2048).
+    output_dir (str | pathlib.Path): Output directory for image files.
+    version (str): Version of resulting image.
+    progress (bool): Check progress bar?
+    **kwargs (dict): Parameters for specific functions.
+
+Returns:
+    SITSCubeModel: A variance data cube.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # create a random forest model
+    rfor_model = sits_train(samples_modis_ndvi, sits_rfor())
+    # create a data cube from local files
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+    # classify a data cube
+    probs_cube = sits_classify(
+        data=cube, ml_model=rfor_model, output_dir=tempfile.gettempdir()
+    )
+    # plot the probability cube
+    plot(probs_cube)
+    # smooth the probability cube using Bayesian statistics
+    var_cube = sits_variance(probs_cube, output_dir=tempfile.gettempdir())
+    # plot the variance cube
+    plot(var_cube)
diff --git a/pysits/docs/content/sits_view.md b/pysits/docs/content/sits_view.md
new file mode 100644
index 0000000..db53524
--- /dev/null
+++ b/pysits/docs/content/sits_view.md
@@ -0,0 +1,145 @@
+View data cubes and samples in leaflet
+
+Uses leaflet to visualize time series, raster cube and classified images.
+
+Args:
+    x (SITSTimeSeriesModel | SITSCubeModel | pandas.DataFrame): time
+        series samples, SOM map, raster cube, probability cube, vector
+        cube, or classified cube.
+    legend (dict): associates labels to colors.
+    palette (str): color palette from RColorBrewer.
+    radius (float): radius of circle markers.
+    add (bool): add image to current leaflet.
+    id_neurons (list[int]): neurons from the SOM map to be shown.
+    band (str): single band for viewing false color images.
+    red (str): band for red color.
+    green (str): band for green color.
+    blue (str): band for blue color.
+    tiles (list[str]): tiles to be plotted (in case of a multi-tile
+        cube).
+    dates (list[str]): dates to be plotted.
+    rev (bool): revert color palette?
+    opacity (float): opacity of segment fill or class cube.
+    max_cog_size (int): maximum size of COG overviews (lines or
+        columns).
+    first_quantile (float): first quantile for stretching images.
+    last_quantile (float): last quantile for stretching images.
+    leaflet_megabytes (float): maximum size for leaflet (in MB).
+    version (str): version name (to compare different classifications).
+    labels (list[str]): labels to be plotted (in case of probs and
+        variance cubes).
+    seg_color (str): color for segment boundaries.
+    line_width (float): line width for segments (in pixels).
+    **kwargs (dict): further specifications for sits_view.
+
+Returns:
+    None: a leaflet object containing either samples or data cubes
+        embedded in a global map that can be visualized directly in a
+        viewer.
+
+Notes:
+    To show a false color image, use "band" to chose one of the bands, "tiles"
+    to select tiles, "first_quantile" and "last_quantile" to set the cutoff
+    points. Choose only one date in the "dates" parameter. The color scheme is
+    defined by either "palette" (use an available color scheme) or legend
+    (user-defined color scheme). To see which palettes are pre-defined, use
+    `cols4all::g4a_gui` or select any ColorBrewer name. The "rev" parameter
+    reverts the order of colors in the palette.
+    To show an RGB composite, select "red", "green" and "blue" bands, "tiles",
+    "dates", "opacity", "first_quantile" and "last_quantile". One can also get
+    an RGB composite, by selecting one band and three dates. In this case, the
+    first date will be shown in red, the second in green and third in blue.
+    Probability cubes are shown in false color. The parameter "labels" controls
+    which labels are shown. If left blank, only the first map is shown. For
+    color control, use "palette", "legend", and "rev" (as described above).
+    Vector cubes have both a vector and a raster component. The vector part are
+    the segments produced by `sits_segment`. Their visual output is controlled
+    by "seg_color" and "line_width" parameters. The raster output works in the
+    same way as the false color and RGB views described above.
+    Classified cubes need information on how to render each class. There are
+    three options: (a) the classes are part of an existing color scheme; (b)
+    the user provides a legend which associates each class to a color; (c) use
+    a generic palette (such as "Spectral") and allocate colors based on this
+    palette. To find out how to create a customized color scheme, read the
+    chapter "Data Visualisation in sits" in the sits book.
+    To compare different classifications, use the "version" parameter to
+    distinguish between the different maps that are shown.
+    Vector classified cubes are displayed as classified cubes, with the
+    segments overlaid on top of the class map, controlled by "seg_color" and
+    "line_width".
+    Samples are shown on the map based on their geographical locations and on
+    the color of their classes assigned in their color scheme. Users can also
+    assign a legend or a palette to choose colors. See information above on the
+    display of classified cubes.
+    For all types of data cubes, the following parameters apply:
+    - opacity: controls the transparency of the map.
+    - max_cog_size: For COG data, controls the level of aggregation to be used
+      for display, measured in pixels, e.g., a value of 512 will select a 512 x
+      512 aggregated image. Small values are faster to show, at a loss of
+      visual quality.
+    - leaflet_megabytes: maximum size of leaflet to be shown associated to the
+      map (in megabytes). Bigger values use more memory.
+    - add: controls whether a new visualisation will be overlaid on top of an
+      existing one. Default is False.
+
+Examples:
+    from pysits import *
+    import tempfile
+
+    # view samples
+    sits_view(cerrado_2classes)
+
+    # create a local data cube
+    data_dir = r_package_dir("extdata/raster/mod13q1", package="sits")
+    modis_cube = sits_cube(
+        source="BDC",
+        collection="MOD13Q1-6.1",
+        data_dir=data_dir
+    )
+
+    # view the data cube
+    sits_view(modis_cube, band="NDVI")
+
+    # train a model
+    rf_model = sits_train(samples_modis_ndvi, sits_rfor())
+
+    # classify the cube
+    modis_probs = sits_classify(
+        data=modis_cube,
+        ml_model=rf_model,
+        output_dir=tempfile.gettempdir()
+    )
+
+    # generate a map
+    modis_label = sits_label_classification(
+        modis_probs,
+        output_dir=tempfile.gettempdir()
+    )
+
+    # view the classified map
+    sits_view(modis_label)
+
+    # add the NDVI band for the first date
+    sits_view(modis_cube,
+        band="NDVI",
+        class_cube=modis_label,
+        dates=sits_timeline(modis_cube)[0],
+        add=True
+    )
+
+    # view the classified map with the RGB image
+    sits_view(modis_cube,
+        red="NDVI", green="NDVI", blue="NDVI",
+        class_cube=modis_label,
+        dates=sits_timeline(modis_cube)[0],
+        add=True
+    )
+
+    # create an uncertainty cube
+    modis_uncert = sits_uncertainty(
+        cube=modis_probs,
+        output_dir=tempfile.gettempdir()
+    )
+
+    # view the uncertainty cube
+    sits_view(modis_uncert, rev=True, add=True)
diff --git a/pysits/docs/content/sits_xgboost.md b/pysits/docs/content/sits_xgboost.md
new file mode 100644
index 0000000..baee0fc
--- /dev/null
+++ b/pysits/docs/content/sits_xgboost.md
@@ -0,0 +1,54 @@
+Train extreme gradient boosting models
+
+This function uses the extreme gradient boosting algorithm. Boosting
+iteratively adds basis functions in a greedy fashion so that each new basis
+function further reduces the selected loss function. This function is a front-
+end to the methods in the "xgboost" package. Please refer to the documentation
+in that package for more details.
+
+Args:
+    samples (SITSTimeSeriesModel): Time series with the training samples.
+    learning_rate (float): Learning rate: scale the contribution of each
+        tree by a factor of 0 < lr < 1 when it is added to the current
+        approximation. Used to prevent overfitting. Default: 0.15
+    min_split_loss (float): Minimum loss reduction to make a further
+        partition of a leaf. Default: 1.
+    max_depth (int): Maximum depth of a tree. Increasing this value makes
+        the model more complex and more likely to overfit. Default: 5.
+    min_child_weight (float): If the leaf node has a minimum sum of
+        instance weights lower than min_child_weight, tree splitting stops.
+        The larger min_child_weight is, the more conservative the algorithm
+        is. Default: 1.
+    max_delta_step (float): Maximum delta step we allow each leaf output to
+        be. If the value is set to 0, there is no constraint. If it is set
+        to a positive value, it can help making the update step more
+        conservative. Default: 1.
+    subsample (float): Percentage of samples supplied to a tree. Default:
+        0.8.
+    nfold (int): Number of the subsamples for the cross-validation.
+    nrounds (int): Number of rounds to iterate the cross-validation
+        (default: 100)
+    nthread (int): Number of threads (default = 6)
+    early_stopping_rounds (int): Training with a validation set will stop
+        if the performance doesn't improve for k rounds.
+    verbose (bool): Print information on statistics during the process
+
+Returns:
+    R: Model fitted to input data (to be passed to `sits_classify`)
+
+Notes:
+    Please refer to the sits documentation available in
+    https://e-sensing.github.io/sitsbook/ for detailed examples.
+
+Examples:
+    from pysits import *
+
+    # Example of training a model for time series classification
+    # Retrieve the samples for Mato Grosso
+    # train a xgboost model
+    ml_model = sits_train(samples_modis_ndvi, ml_method=sits_xgboost)
+    # classify the point
+    point_ndvi = sits_select(point_mt_6bands, bands="NDVI")
+    # classify the point
+    point_class = sits_classify(data=point_ndvi, ml_model=ml_model)
+    plot(point_class)
diff --git a/pysits/docs/content/summary.md b/pysits/docs/content/summary.md
new file mode 100644
index 0000000..64ba928
--- /dev/null
+++ b/pysits/docs/content/summary.md
@@ -0,0 +1,36 @@
+Summarize sits objects and data cubes.
+
+This is a generic function that produces a summary of the input object.
+The behavior and available parameters depend on the specific type of
+object provided. It can summarize time series, raster cubes, classified
+cubes, variance cubes, and accuracy assessments.
+
+Args:
+    object (SITSTimeSeriesModel | SITSCubeModel | SITSConfusionMatrix):
+            The object to be summarized. Supported types include: a set
+            of time series; a raster data cube; a classified data cube;
+            a variance data cube; an accuracy matrix for training data;
+            or an accuracy matrix for area data.
+    **kwargs (dict): Further specifications for the summary. The
+            accepted keyword arguments depend on the type of `object`:
+
+            For a raster data cube:
+                tile: Tile to be summarized.
+                date: Date to be summarized.
+
+            For a variance data cube:
+                intervals: Intervals to calculate the quantiles.
+                sample_size: The approximate size of samples that
+                    will be extracted from the variance cube (by
+                    tile).
+                multicores: Number of cores to summarize data
+                    (min = 1, max = 2048).
+                memsize: Memory in GB available to summarize data
+                    (min = 1, max = 16384).
+                quantiles: Quantiles to be shown.
+
+Returns:
+    str: A summary of the input object. Depending on the input type,
+        this is a summary of the time series, the raster data cube, the
+        classified data cube, the variance data cube, or the
+        sample/area accuracy.