Module analysis

All functions here are automatically loaded with from sits import analysis.

analysis.initialize_dask_client

sits.analysis.initialize_dask_client(n_cores=False, threads_per_worker=1)

Initializes a Dask distributed client using by default all available CPU cores, with one thread per worker. Restarts the client to ensure a clean state.

Parameters:

  • n_cores (bool, optional) – Number of CPU cores. Defaults to False (i.e. all available CPU cores).

  • crs_out (int, optional) – Number of threads per worker. Defaults to 1.

Returns:

A Dask distributed client instance.

Return type:

Client

Example

>>> client = initialize_dask_client()

analysis.date_range

sits.analysis.date_range(start_date, end_date, freq='D')

Generate a range of dates between two given dates with a specified frequency.

Parameters:

  • start_date (str) – start date ‘yyyy-mm-dd’ of the range.

  • end_date (str) – end date ‘yyyy-mm-dd’ of the range.

  • freq (str, optional) – accepts any valid pandas frequency string. Default to ‘D’ for daily.

Returns:

pandas DatetimeIndex object representing the date range.

Return type:

date_range (pd.DatetimeIndex)

Examples

>>> ts_predict = date_range('2025-01-01', '2025-07-01', freq='D')

analysis.reindexTS

sits.analysis.reindexTS(df, freq='D', regular_freq=False, interpolate=False, method='linear')

Reindex a time series DataFrame to a regular frequency, with optional aggregation and interpolation.

Parameters:

  • df (pd.DataFrame) – input DataFrame with a datetime index.

  • freq (str, optional) – desired frequency for the output time series. Default to ‘D’ for daily.

  • regular_freq (bool, optional) – If True, enforces regular frequency. If freq is not ‘D’, this must be False, otherwise a ValueError is raised. Defaults to False.

  • interpolate (bool, optional) – If True, interpolates missing values using the specified method. Defaults to False.

  • method (str, optional) – interpolation method to use. Accepts any method supported by pandas.DataFrame.interpolate. Defaults to ‘linear’.

Returns:

a DataFrame reindexed to the specified frequency,

with optional interpolation applied.

Return type:

df (pd.DataFrame)

Raises:

ValueError – if regular_freq is True while freq is not ‘D’.

Notes

  • The function uses internal helpers: _ensure_datetime_index, _resample_df, _convert_to_period, _regularize_index, and _fill_nan.

  • The index is normalized to remove time components before processing.

Examples

>>> ts_train = reindexTS(df,
...     freq='D',
...     regular_freq=True,
...     interpolate=True,
...     method='linear')

analysis.sktime_fitpred

sits.analysis.sktime_fitpred(ts, time_index, predict_time, model_name='ThetaForecaster', model_params=None, freq='D', reindex_params=None)

Fit a time series forecasting model from sktime and predict future values.

This function prepares a time series dataset, optionally reindexes it to a regular frequency, fits a forecasting model from the sktime library, and returns predictions for specified future timestamps.

Parameters:

  • ts (array-like) – the time series values to be modeled.

  • time_index (array-like) – the corresponding datetime index for the time series values.

  • predict_time (array-like) – a list or array of future timestamps for which predictions are required.

  • model_name (str, optional) – name of the sktime forecasting model to use. Must be a valid forecaster registered in sktime. Defaults to ‘ThetaForecaster’.

  • model_params (dict, optional) – dictionary of parameters to initialize the forecasting model. If None, default parameters are used.

  • freq (str, optional) – frequency string used to regularize the time index. Default to ‘D’ for daily.

  • reindex_params (dict, optional) – additional keyword arguments passed to the reindexTS function for time series reindexing.

Returns:

array of predicted values corresponding to predict_time.

Return type:

np.ndarray

Raises:

  • ValueError if the specified model_name is not found in –

  • the sktime registry. –

Examples

>>> predict = sktime_fitpred(arr_train, time_train, time_pred)

analysis.xr_forecast

sits.analysis.xr_forecast(dataarray, predict_time, model_name='ThetaForecaster', model_params=None, reindex_params=None, freq='D')

Apply a sktime forecasting model to each pixel of a time series DataArray.

This function performs pixel-wise time series forecasting on an xarray.DataArray using a specified sktime model. It applies the sktime_fitpred function across the spatial dimensions using xarray.apply_ufunc, enabling parallelized computation with Dask.

Parameters:

  • dataarray (xr.DataArray) – a 3D xarray DataArray with dimensions (‘time’, ‘y’, ‘x’), representing the time series data for each pixel.

  • predict_time (array-like) – a list or array of future timestamps for which predictions are required.

  • model_name (str, optional) – name of the sktime forecasting model to use. Must be a valid forecaster registered in sktime. Default to ‘ThetaForecaster’.

  • model_params (dict, optional) – dictionary of parameters to initialize the forecasting model. If None, default parameters are used.

  • reindex_params (dict, optional) – additional keyword arguments passed to the reindexTS function for time series reindexing.

  • freq (str, optional) – frequency string used to regularize the time index. Default to ‘D’ for daily.

Returns:

a DataArray with dimensions (‘time’, ‘y’, ‘x’) containing the forecasted values for each pixel.

Return type:

result (xr.DataArray)

Notes

  • The function uses Dask for parallelized execution, making it suitable

for large datasets. - The output time coordinate is replaced with predict_time.

Examples

>>> result = xr_forecast(dataarray, predict_time)

analysis.ClearCut

class sits.analysis.ClearCut(da: DataArray)

Bases: object

This class aims to detect changes in univariate time series (e.g., spectral indices such as NDVI). It identifies potential change points by comparing values in two moving windows:

  • A backward window summarizes past observations.

  • A forward window summarizes upcoming observations.

A change is flagged when the forward window shows a significant drop relative to the backward window, according to a configurable threshold. This method is an adaptation of the OTB ClearCutsMultitemporalDetection application (R. Cresson, INRAE).

Typical use case: detecting vegetation loss or disturbance events in spectral index time series (NDVI, EVI, etc.), where sudden drops may indicate deforestation, fire, or degradation.

da

index time series (‘time’, ‘band’, ‘y’, ‘x’)

Type:

xr.Dataarray

Parameters:

da (xr.Dataarray) – index time series (‘time’, ‘band’, ‘y’, ‘x’)

Example

>>> ndvi_ts = ClearCut(Dataarray)
detect_anomalies(thresholds: list = [0.2, 0.3, 0.4], anomaly_type: str = 'absolute', window_backward: int = 720, window_forward: int = 60, min_obs_backward: int = 5, min_obs_forward: int = 2, out_crs: str = 'epsg:3035', store_magnitude: bool = False)

Detect anomalies in a univariate time series (e.g., spectral indices such as NDVI) by comparing values between backward and forward moving windows.

For each date in the series (excluding margins defined by the window sizes), the method computes:

  • Mean of the backward window (past observations).

  • Mean of the forward window (future observations).

  • Magnitude of change as the absolute difference between the two means.

An anomaly is flagged when the magnitude exceeds the first threshold. The function records the first, maximum, and last anomalies detected, along with their dates, and classifies anomalies based on the maximum magnitude relative to the provided thresholds.

Parameters:

  • thresholds (list, optional) – list of float. Thresholds used to classify anomaly magnitudes. Defaults to [0.2, 0.3, 0.4].

  • anomaly_type (str) – detection direction. Defaults to “absolute”. - “absolute”: detects any change > threshold (default). - “drop”: detects only negative changes (e.g., vegetation loss). - “increase”: detects only positive changes (e.g., regrowth).

  • window_backward (int, optional) – size of the backward window in days. Defaults to 720.

  • window_forward (int, optional) – size of the forward window in days. Defaults to 60.

  • min_obs_backward (int, optional) – minimum number of valid observations required in the backward window. Defaults to 5.

  • min_obs_forward (int, optional) – minimum number of valid observations required in the forward window. Defaults to 2.

  • out_crs (str, optional) – Coordinate Reference System (CRS) to assign to the output dataset. Defaults to “epsg:3035”.

  • store_magnitude (bool, optional) – Export of the magnitude timeseries. Defaults to False.

Returns:

Results are stored with the following variables:

  • magnitude_first (float): Magnitude of the first anomaly detected.

  • date_first (Unix epoch): Date (days since epoch) of the first anomaly.

  • magnitude_max (float): Maximum anomaly magnitude detected.

  • date_max (Unix epoch): Date of the maximum anomaly.

  • magnitude_last (float): Magnitude of the last anomaly detected.

  • date_last (Unix epoch): Date of the last anomaly.

  • mask (binary): Boolean mask indicating anomaly presence.

  • inrange (binary): Boolean mask indicating whether sufficient observations were available in both windows.

  • classif (int): Classification layer based on thresholds.

ClearCut.magnitude_ts (xr.Dataarray, optional): datacube of magnitudes

Return type:

ClearCut.detection (xr.Dataset)

Notes

  • Dates are stored as integer days since 1970-01-01.

  • The classification is performed using the maximum anomaly magnitude compared against the provided thresholds.

  • The output dataset is tagged with the specified CRS for spatial consistency in geospatial workflows.

Example

>>> ndvi_ts.detect_anomalies()
season_of_interest(start: str = '06-01', end: str = '08-31', compute=True)

Filters the datacube to retain only dates between start and end (format: ‘mm-dd’).

Parameters:

  • start (str, optional) – start date in ‘mm-dd’ format. Defaults to ‘06-01’.

  • end (str, optional) – end date in ‘mm-dd’ format. Defaults to ‘08-31’.

Returns:

Filtered Dataarray with only dates in the specified range.

Return type:

ClearCut.da (xr.Dataarray)

Example

>>> ndvi_ts.season_of_interest()

analysis.sieve_maj

sits.analysis.sieve_maj(dataarray, min_size=3, window_size=3, ignore_nan=True, connectivity=1, out_crs='epsg:3035')

Remove small connected objects and replace them with local majority value. Preserves original coordinates and CRS.

Parameters:

  • dataarray (xr.DataArray) – input array with integer classification values.

  • min_size (int, optional) – minimum size (number of pixels) to keep. Defaults to 3.

  • window_size (int, optional) – size of the moving window for majority filter. Defaults to 3

  • ignore_nan (bool, optional) – if True, NaNs are treated as background (converted to 0). Defaults to True.

  • connectivity (int, optional) – connectivity for labeling (1=4-connectivity, 2=8-connectivity). Defaults to 1

  • out_crs (str, optional) – output CRS. Defaults to ‘epsg:3035’.

Returns:

filtered DataArray with original coords and CRS.

Return type:

xr.DataArray

Example

>>> sieve_maj(ndvi_ts.detection.classif)

analysis.SitsPlotter

class sits.analysis.SitsPlotter(sits_object, band_index, daterange, break_date=datetime.datetime(2024, 1, 1, 0, 0), i=0, j=0, coords=None, ylim=[-1, 1], filter_method=None, filter_value=1.5)

Bases: object

Visualization utility for Satellite Image Time Series (SITS) data.

This class handles extraction, outlier filtering, and plotting of temporal data for a specific spatial coordinate. It segments the time series into ‘fitting’ and ‘monitoring’ windows relative to a break date and supports statistical outlier removal (Sigma or IQR).

Key features: * Time Series Segmentation: Splits data into training (fit) and

monitoring windows around a specified break date.

  • Statistical Filtering: Supports Sigma (mean/std) or IQR-based outlier removal.

  • Detection Integration: Interfaces with change detection algorithms (e.g., ClearCut) to visualize temporal anomalies and magnitudes.

data

The extracted time series for the specified coordinate.

Type:

xr.DataArray

Parameters:

  • sits_object (xr.Dataset) – Time series dataset (e.g., StacAttack.cube or StacAttack.indices).

  • band_index (str) – The specific band or index name to display (e.g., ‘B04’, ‘NDVI’).

  • daterange (list[datetime]) – Period to display as [start_date, end_date].

  • break_date (datetime, optional) – The transition point between fitting and monitoring periods. Defaults to 2024-01-01.

  • i (int, optional) – Image column index (x-axis in pixel space). Defaults to 0.

  • j (int, optional) – Image row index (y-axis in pixel space). Defaults to 0.

  • coords (list[float], optional) – Spatial coordinates [x, y] in the same CRS as sits_object. If provided, overrides i and j. Defaults to None.

  • ylim (list[float], optional) – Y-axis limits for the plot. Defaults to [-1, 1].

  • filter_method (str, optional) – Outlier detection method: ‘sigma’ or ‘iqr’. Defaults to None.

  • filter_value (float, optional) – Threshold for filtering (number of std devs for ‘sigma’ or multiplier for ‘iqr’). Defaults to 1.5.

Example

>>> from datetime import datetime
>>> daterange = [datetime(2023, 1, 1), datetime(2025, 1, 1)]
>>> plotter = SitsPlotter(
...     sits_object=stacObj.indices,
...     band_index='NDMI',
...     daterange=daterange,
...     break_date=datetime(2024, 1, 1),
...     filter_method='sigma'
... )
>>> plotter.plot_model(degree=3,
...     include_trend=True,
...     envelope_quantiles=(0.05, 0.95)
... )
>>> plotter.run_detection(thresholds=[0.2, 0.3, 0.4])
>>> plotter.save_plot("sits_out/ts_plot.png")
plot_model(degree=1, include_trend=True, color='crimson', envelope_quantiles=(0.25, 0.75))

Fits a regression model to the fitting window and plots the projection.

This method uses Ordinary Least Squares (OLS) to fit a model to the ‘fit_window’ data. The model can include a linear trend and harmonic components (via the design matrix). It then extrapolates this model across the full ‘daterange’ to visualize the expected behavior versus actual observations.

The uncertainty envelope is calculated based on the distribution of residuals within the training (fitting) period.

Parameters:

  • degree (int, optional) – The polynomial degree or number of harmonic components for the model. Defaults to 1.

  • include_trend (bool, optional) – Whether to include a linear trend component in the model. Defaults to True.

  • color (str, optional) – Color of the model line and uncertainty envelope. Defaults to ‘crimson’.

  • envelope_quantiles (tuple[float, float], optional) – The lower and upper quantiles of the residuals used to define the shaded uncertainty area. Defaults to (0.25, 0.75).

Returns:

The fitted model coefficients (beta weights).

Return type:

SitsPlotter.beta (np.ndarray)

Example

>>> plotter.plot_model(degree=3,
...     include_trend=True,
...     envelope_quantiles=(0.05, 0.95)
... )
run_detection(detector_class=<class 'sits.analysis.ClearCut'>, thresholds=[0.2, 0.3, 0.4], anomaly_type='absolute', **kwargs)

Runs change detection on the time series and visualizes the results.

This method initializes a detector (e.g., ClearCut), processes the time series for anomalies, and dynamically appends a magnitude subplot to the main figure if one does not already exist. It identifies the maximum anomaly date and marks it across both the main index plot and the magnitude plot.

Parameters:

  • detector_class (class, optional) – The class to use for anomaly detection. Must implement .detect_anomalies() and store magnitude results. Defaults to ClearCut.

  • thresholds (list[float], optional) – Values used to identify significant changes. All values must have the same sign (either all positive or all negative). Defaults to [0.2, 0.3, 0.4].

  • anomaly_type (str, optional) – The mode for anomaly detection. Defaults to ‘absolute’. - “absolute”: detects any change > threshold (default). - “drop”: detects only negative changes. - “increase”: detects only positive changes.

  • **kwargs – Additional arguments passed directly to the detector’s detect_anomalies method.

Example

>>> plotter.run_detection(thresholds=[0.2, 0.3, 0.4])
save_plot(output_path='output', filename=None)

Saves the current figure to a specified file path.

Automatically creates the output directory if it does not exist. Uses a high-resolution DPI and includes the external legend in the final file.

Parameters:

  • output_path (str, optional) – The directory where the plot will be saved. Defaults to “output”.

  • filename (str, optional) – The name of the file (e.g., “plot.png”). If None, defaults to “sits_{band_index}.png”.

Returns:

The full absolute path to the saved image file.

Return type:

str