Merge pull request #1127 from M3nin0/feature/detections-api

detect changes: introduce detect changes api

DESCRIPTION

@@ -130,7 +130,9 @@ Collate:
  'api_cube.R'
  'api_data.R'
  'api_debug.R'
+ 'api_detect_changes.R'
  'api_download.R'
+ 'api_dtw.R'
  'api_environment.R'
  'api_factory.R'
  'api_file_info.R'
@@ -215,6 +217,9 @@ Collate:
  'sits_cube_copy.R'
  'sits_clean.R'
  'sits_cluster.R'
+ 'sits_detect_change.R'
+ 'sits_detect_change_method.R'
+ 'sits_dtw.R'
  'sits_factory.R'
  'sits_filters.R'
  'sits_geo_dist.R'

NAMESPACE

@@ -338,6 +338,8 @@ S3method(sits_cube,default)
 S3method(sits_cube,local_cube)
 S3method(sits_cube,sar_cube)
 S3method(sits_cube,stac_cube)
+S3method(sits_detect_change,default)
+S3method(sits_detect_change,sits)
 S3method(sits_get_data,csv)
 S3method(sits_get_data,data.frame)
 S3method(sits_get_data,default)
@@ -441,6 +443,9 @@ export(sits_config)
 export(sits_config_show)
 export(sits_cube)
 export(sits_cube_copy)
+export(sits_detect_change)
+export(sits_detect_change_method)
+export(sits_dtw)
 export(sits_factory_function)
 export(sits_filter)
 export(sits_formula_linear)

R/RcppExports.R

@@ -9,6 +9,10 @@ weighted_uncert_probs <- function(data_lst, unc_lst) {
  .Call(`_sits_weighted_uncert_probs`, data_lst, unc_lst)
 }
 
+dtw_distance <- function(ts1, ts2) {
+ .Call(`_sits_dtw_distance`, ts1, ts2)
+}
+
 C_kernel_median <- function(x, ncols, nrows, band, window_size) {
  .Call(`_sits_C_kernel_median`, x, ncols, nrows, band, window_size)
 }

R/api_detect_changes.R

@@ -0,0 +1,44 @@
+#' @title Detect changes in time-series using various methods.
+#' @name .detect_change_ts
+#' @keywords internal
+#' @noRd
+.detect_change_ts <- function(samples,
+ cd_method,
+ filter_fn,
+ multicores,
+ progress) {
+ # Start parallel workers
+ .parallel_start(workers = multicores)
+ on.exit(.parallel_stop(), add = TRUE)
+ # Get bands from model
+ bands <- .ml_bands(cd_method)
+ # Update samples bands order
+ if (any(bands != .samples_bands(samples))) {
+ samples <- .samples_select_bands(samples = samples,
+ bands = bands)
+ }
+ # Apply time series filter
+ if (.has(filter_fn)) {
+ samples <- .apply_across(data = samples,
+ fn = filter_fn)
+ }
+ # Divide samples in chunks to parallel processing
+ parts <- .pred_create_partition(pred = samples, partitions = multicores)
+ # Detect changes!
+ detections <- .jobs_map_parallel_dfr(parts, function(part) {
+ # Get samples
+ values <- .pred_part(part)
+ # Detect changes! For detection, models can be time-aware. So, the
+ # complete data, including dates, must be passed as argument.
+ detections <- cd_method(values[["time_series"]])
+ detections <- tibble::tibble(detections)
+ # Prepare result
+ result <- tibble::tibble(data.frame(values, detections = detections))
+ class(result) <- class(values)
+ # return
+ result
+
+ }, progress = progress)
+
+ return(detections)
+}

R/api_dtw.R

@@ -0,0 +1,56 @@
+
+#' @title Extract temporal pattern of samples using temporal median.
+#' @name .pattern_temporal_median
+#' @keywords internal
+#' @noRd
+.pattern_temporal_median <- function(samples) {
+ samples |>
+ dplyr::group_by(.data[["label"]]) |>
+ dplyr::group_map(function(data, name) {
+ ts_median <- dplyr::bind_rows(data[["time_series"]]) |>
+ dplyr::group_by(.data[["Index"]]) |>
+ dplyr::summarize(dplyr::across(dplyr::everything(),
+ stats::median, na.rm = TRUE)) |>
+ dplyr::select(-.data[["Index"]])
+
+ ts_median["label"] <- name
+ ts_median
+ })
+}
+
+#' @title Calculate the DTW distance between label patterns and sample data.
+#' @name .pattern_distance_dtw
+#' @description This function calculates the DTW distance between label patterns
+#' and sample data in a given temporal window.
+#' @keywords internal
+#' @noRd
+.pattern_distance_dtw <- function(data, patterns, windows) {
+ # Calculate the DTW distance between `data` and `patterns`
+ purrr::map_dfc(1:length(patterns), function(pattern_index) {
+ # Get pattern metadata
+ pattern <- patterns[pattern_index][[1]]
+ pattern_label <- unique(pattern[["label"]])
+ # Get pattern data
+ pattern_ts <- dplyr::select(pattern, -.data[["label"]])
+ pattern_ts <- as.matrix(pattern_ts)
+ # Windowed search
+ distances <- purrr::map_df(windows, function(window) {
+ # Get time-series in the window
+ data_in_window <-
+ dplyr::filter(data,
+ .data[["Index"]] >= window[["start"]],
+ .data[["Index"]] <= window[["end"]])
+ # Remove the time reference column
+ data_in_window <- dplyr::select(data_in_window, -.data[["Index"]])
+ # Transform values in matrix (as expected in the cpp code)
+ data_in_window <- as.matrix(data_in_window)
+ data_in_window <- data_in_window
+ # Calculate distance
+ distance_from_pattern <- dtw_distance(data_in_window, pattern_ts)
+ # Prepare result and return it
+ data.frame(distance = distance_from_pattern)
+ })
+ # Associate the pattern name with the distances
+ stats::setNames(distances, pattern_label)
+ })
+}

R/api_period.R

@@ -47,3 +47,37 @@ NULL
  unit <- c(D = "day", M = "month", Y = "year")
  unit[[gsub("^P[0-9]+([DMY])$", "\\1", period)]]
 }
+
+#' @describeIn period_api Create period windows.
+#' @returns \code{.period_windows()}: Period windows.
+#' @noRd
+.period_windows <- function(period, step, start_date, end_date) {
+ # Transform `period` and `step` strings in duration
+ period_duration <- lubridate::as.duration(period)
+ step_duration <- lubridate::as.duration(step)
+ # Transform `start_date` and `end_date` to date
+ start_date <- as.Date(start_date)
+ end_date <- as.Date(end_date)
+ # Final period windows
+ period_windows <- list()
+ # Define first time period (used as part of the step)
+ current_start <- start_date
+ # Create period windows
+ while(current_start < end_date) {
+ # Create the window: current start date + step
+ current_end <- current_start + period_duration
+ # Avoid window definition beyond the end date
+ if (current_end > end_date) {
+ current_end <- end_date
+ }
+ # Save period window
+ period_windows <-
+ c(period_windows, list(c(
+ start = as.Date(current_start),
+ end = as.Date(current_end)
+ )))
+ # Move to the next window date: current start date + step
+ current_start <- current_start + step_duration
+ }
+ period_windows
+}

R/sits_detect_change.R

@@ -0,0 +1,60 @@
+#' @title Detect changes in time series
+#' @name sits_detect_change
+#'
+#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
+#' @author Felipe Carlos, \email{efelipecarlos@@gmail.com}
+#' @author Felipe Carvalho, \email{felipe.carvalho@@inpe.br}
+#'
+#' @description Given a set of time series or an image, this function compares
+#' each time series with a set of change/no-change patterns, and indicates
+#' places and dates where change has been detected.
+#'
+#' @param data Set of time series.
+#' @param cd_method Change detection method (with parameters).
+#' @param ... Other relevant parameters.
+#' @param filter_fn Smoothing filter function to be applied to the data.
+#' @param multicores Number of threads to process the time series.
+#' @param progress Show progress bar?
+#' @return Set of time series where significant change has been
+#' detected.
+#' @export
+sits_detect_change <- function(data,
+ cd_method,
+ ...,
+ filter_fn = NULL,
+ multicores = 2L,
+ progress = TRUE) {
+ UseMethod("sits_detect_change", data)
+}
+
+#' @rdname sits_detect_change
+#' @export
+sits_detect_change.sits <- function(data,
+ cd_method,
+ ...,
+ filter_fn = NULL,
+ multicores = 2L,
+ progress = TRUE) {
+ # set caller for error messages
+ .check_set_caller("sits_detect_change_sits")
+ # Pre-conditions
+ data <- .check_samples_ts(data)
+ .check_is_sits_model(cd_method)
+ .check_int_parameter(multicores, min = 1, max = 2048)
+ .check_progress(progress)
+ # Do detection
+ detections <- .detect_change_ts(
+ samples = data,
+ cd_method = cd_method,
+ filter_fn = filter_fn,
+ multicores = multicores,
+ progress = progress
+ )
+ return(detections)
+}
+
+#' @rdname sits_detect_change
+#' @export
+sits_detect_change.default <- function(data, cd_method, ...) {
+ stop("Input should be a sits tibble or a data cube")
+}

R/sits_detect_change_method.R

@@ -0,0 +1,35 @@
+#' @title Create detect change method.
+#' @name sits_detect_change_method
+#'
+#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
+#' @author Felipe Carlos, \email{efelipecarlos@@gmail.com}
+#'
+#' @description Prepare detection change method. Currently, sits supports the
+#' following methods: 'dtw' (see \code{\link[sits]{sits_dtw}})
+#'
+#' @param samples Time series with the training samples.
+#' @param cd_method Change detection method.
+#' @return Change detection method prepared
+#' to be passed to
+#' \code{\link[sits]{sits_detect_change}}
+#' @export
+#'
+sits_detect_change_method <- function(samples, cd_method = sits_dtw()) {
+ # set caller to show in errors
+ .check_set_caller("sits_detect_change_method")
+ # check if samples are valid
+ .check_samples_train(samples)
+ # is the train method a function?
+ .check_that(inherits(cd_method, "function"),
+ msg = .conf("messages", "sits_detect_change_method_model")
+ )
+ # are the timelines OK?
+ timeline_ok <- .timeline_check(samples)
+ .check_that(timeline_ok,
+ msg = .conf("messages", "sits_detect_change_method_timeline")
+ )
+ # compute the training method by the given data
+ result <- cd_method(samples)
+ # return a valid detect change method
+ return(result)
+}

R/sits_dtw.R

@@ -0,0 +1,111 @@
+#' @title Dynamic Time Warping for Detect changes.
+#' @name sits_dtw
+#'
+#' @author Felipe Carlos, \email{efelipecarlos@@gmail.com}
+#' @author Felipe Carvalho, \email{felipe.carvalho@@inpe.br}
+#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
+#' @author Rolf Simoes, \email{rolf.simoes@@inpe.br}
+#'
+#' @description Create a Dynamic Time Warping (DTW) method for the
+#' \code{\link[sits]{sits_detect_change_method}}.
+#'
+#' @param samples Time series with the training samples.
+#' @param ... Other relevant parameters.
+#' @param threshold Threshold used to define if an event was detected.
+#' Default is `Inf`.
+#' @param window ISO8601-compliant time period used to define the
+#' DTW moving window, with number and unit,
+#' where "D", "M" and "Y" stands for days, month and
+#' year; e.g., "P16D" for 16 days.
+#' @return Change detection method prepared to be passed to
+#' \code{\link[sits]{sits_detect_change_method}}
+#' @export
+#'
+sits_dtw <-
+ function(samples = NULL,
+ ...,
+ threshold = Inf,
+ window = NULL) {
+ .check_set_caller("sits_dtw")
+ train_fun <-
+ function(samples) {
+ # Check parameters
+ .check_period(window)
+ .check_null_parameter(threshold)
+ # Sample labels
+ labels <- .samples_labels(samples)
+ # Get samples patterns (temporal median)
+ train_samples_patterns <- .pattern_temporal_median(samples)
+ # Define detection function
+ detect_change_fun <- function(values) {
+ # Extract dates
+ dates <- values[[1]][["Index"]]
+ dates_min <- min(dates)
+ dates_max <- max(dates)
+ # Assume time-series are regularized, then use the period
+ # as the step of the moving window. As a result, we have
+ # one step per iteration.
+ dates_step <- lubridate::as.period(
+ lubridate::int_diff(dates)
+ )
+ dates_step <- dates_step[[1]]
+ # Create comparison windows
+ comparison_windows <- .period_windows(
+ period = window,
+ step = dates_step,
+ start_date = dates_min,
+ end_date = dates_max
+ )
+ # Do the change detection for each time-series
+ purrr::map(values, function(value_row) {
+ # Search for the patterns
+ patterns_distances <- .pattern_distance_dtw(
+ data = value_row,
+ patterns = train_samples_patterns,
+ windows = comparison_windows
+ )
+ # Remove distances out the user-defined threshold
+ patterns_distances[patterns_distances > threshold] <- NA
+ # Define where each label was detected. For this, first
+ # get from each label the minimal distance
+ detections_idx <-
+ apply(patterns_distances, 2, which.min)
+ detections_name <- names(detections_idx)
+ # For each label, extract the metadata where they had
+ # minimal distance
+ purrr::map_df(1:length(detections_idx), function(idx) {
+ # Extract detection name and inced
+ detection_name <- detections_name[idx]
+ detection_idx <- detections_idx[idx]
+ # Extract detection distance (min one defined above)
+ detection_distance <-
+ patterns_distances[detection_idx,]
+ detection_distance <-
+ detection_distance[detection_name]
+ detection_distance <-
+ as.numeric(detection_distance)
+ # Extract detection dates
+ detection_dates <-
+ comparison_windows[[detection_idx]]
+ # Prepare result and return it!
+ data.frame(
+ change = detection_name,
+ distance = detection_distance,
+ from = detection_dates[["start"]],
+ to = detection_dates[["end"]]
+ )
+ })
+ })
+ }
+ # Set model class
+ detect_change_fun <- .set_class(detect_change_fun,
+ "dtw_model",
+ "sits_model",
+ class(detect_change_fun))
+ return(detect_change_fun)
+ }
+ # If samples is informed, train a model and return a predict function
+ # Otherwise give back a train function to train model further
+ result <- .factory_function(samples, train_fun)
+ return(result)
+ }