diff --git a/pysteps/nowcasts/anvil.py b/pysteps/nowcasts/anvil.py
index 9da0fb47e..c9dfbe774 100644
--- a/pysteps/nowcasts/anvil.py
+++ b/pysteps/nowcasts/anvil.py
@@ -23,6 +23,7 @@
 from scipy.ndimage import gaussian_filter
 from pysteps import cascade, extrapolation
 from pysteps.nowcasts.utils import nowcast_main_loop
+from pysteps.nowcasts import utils as nowcast_utils
 from pysteps.timeseries import autoregression
 from pysteps import utils
 
@@ -184,7 +185,7 @@ def forecast(
         starttime_init = time.time()
 
     m, n = vil.shape[1:]
-    vil = vil.copy()
+    vil = vil[-(ar_order+2) :].copy()
 
     if rainrate is None and apply_rainrate_mask:
         rainrate_mask = vil[-1, :] < 0.1
@@ -273,6 +274,10 @@ def worker(vil, i):
                 gamma[i, 0, :], gamma[i, 1, :]
             )
 
+    nowcast_utils.print_corrcoefs(
+        np.nanmean(gamma, axis=(-2,-1))
+        )
+    
     # estimate the parameters of the ARI models
     phi = []
     for i in range(n_cascade_levels):
@@ -284,6 +289,10 @@ def worker(vil, i):
             phi_ = _estimate_ar1_params(gamma[i, :])
         phi.append(phi_)
 
+    nowcast_utils.print_ar_params(
+        np.nanmean(np.array(phi), axis=(-2,-1))
+        )
+    
     vil_dec = vil_dec[:, -(ar_order + 1) :, :]
 
     if measure_time:
@@ -339,7 +348,7 @@ def _check_inputs(vil, rainrate, velocity, timesteps, ar_order):
                 "rainrate.shape = %s, but a two-dimensional array expected"
                 % str(rainrate.shape)
             )
-    if vil.shape[0] != ar_order + 2:
+    if vil.shape[0] < ar_order + 2:
         raise ValueError(
             "vil.shape[0] = %d, but vil.shape[0] = ar_order + 2 = %d required"
             % (vil.shape[0], ar_order + 2)
diff --git a/pysteps/verification/detcatscores.py b/pysteps/verification/detcatscores.py
index 83a8d1512..6162fbeed 100644
--- a/pysteps/verification/detcatscores.py
+++ b/pysteps/verification/detcatscores.py
@@ -18,9 +18,9 @@
 
 import collections
 import numpy as np
+from pysteps.verification.spatialscores import smooth_fields, normalize_forecast_obs
 
-
-def det_cat_fct(pred, obs, thr, scores="", axis=None):
+def det_cat_fct(pred, obs, thr, scores="", axis=None, return_table=False):
     """
     Calculate simple and skill scores for deterministic categorical
     (dichotomous) forecasts.
@@ -94,7 +94,7 @@ def det_cat_fct(pred, obs, thr, scores="", axis=None):
 
     contab = det_cat_fct_init(thr, axis)
     det_cat_fct_accum(contab, pred, obs)
-    return det_cat_fct_compute(contab, scores)
+    return det_cat_fct_compute(contab, scores, return_table)
 
 
 def det_cat_fct_init(thr, axis=None):
@@ -198,11 +198,22 @@ def det_cat_fct_accum(contab, pred, obs):
         obs = obs[None, :]
         axis = (0,)
     axis = tuple([a for a in axis if a >= 0])
-
+    
+    # Auto-detect binary inputs: unique values subset of {0,1} or {False,True}
+    def is_binary(arr):
+        vals = np.unique(arr[~np.isnan(arr)])
+        return np.all(np.isin(vals, [0, 1, False, True]))
+    
+    binary_input = is_binary(pred) and is_binary(obs)
+    
     # apply threshold
-    predb = pred > contab["thr"]
-    obsb = obs > contab["thr"]
-
+    if binary_input:
+        predb = pred.astype(bool)
+        obsb = obs.astype(bool)
+    else:
+        predb = pred > contab["thr"]
+        obsb = obs > contab["thr"]
+    
     # calculate hits, misses, false positives, correct rejects
     H_idx = np.logical_and(predb == 1, obsb == 1)
     F_idx = np.logical_and(predb == 1, obsb == 0)
@@ -263,7 +274,7 @@ def det_cat_fct_merge(contab_1, contab_2):
     return contab
 
 
-def det_cat_fct_compute(contab, scores=""):
+def det_cat_fct_compute(contab, scores="", return_table=False):
     """
     Compute simple and skill scores for deterministic categorical
     (dichotomous) forecasts from a contingency table object.
@@ -332,8 +343,14 @@ def get_iterable(x):
     M = 1.0 * contab["misses"]  # false negatives
     F = 1.0 * contab["false_alarms"]  # false positives
     R = 1.0 * contab["correct_negatives"]  # true negatives
-
+    
     result = {}
+    if return_table:
+        result["H"] = H
+        result["M"] = M
+        result["F"] = F
+        result["R"] = R
+        
     for score in scores:
         # catch None passed as score
         if score is None:
@@ -401,3 +418,78 @@ def get_iterable(x):
             result["F1"] = F1
 
     return result
+
+
+def det_cat_fct_spatial(pred, obs, thr, scales, scores="", axis=None, return_table=False, smooth_thr=0.01):
+    """
+    Perform neighborhood-based categorical verification across thresholds and scales.
+
+    Parameters
+    ----------
+    pred : np.ndarray
+        Forecast array. Shape can be:
+        - (n_members, n_leadtimes, y, x)
+        - (n_leadtimes, y, x)
+        - (n_members, y, x)
+        - (y, x)
+    obs : np.ndarray
+        Observation array. Shape can be:
+        - (n_leadtimes, y, x)
+        - (y, x)
+    thr : float or list of float
+        Threshold(s) for binarization.
+    scales : list of int
+        Smoothing window sizes.
+    axis : int or tuple of int, optional
+        Axes over which to aggregate contingency table.
+    scores : str or list of str, optional
+        Score names to compute. If empty, all supported scores are used.
+    smooth_thr : float
+        Threshold to re-binarize smoothed fields.
+
+    Returns
+    -------
+    results : dict
+        Dictionary of score arrays:
+        results[score_name] = np.ndarray of shape (n_thresholds, n_scales, n_members, n_leadtimes)
+    """
+    thr = np.atleast_1d(thr)
+    scales = np.atleast_1d(scales)
+
+    # Normalize forecast and observation shapes
+    pred, obs = normalize_forecast_obs(pred, obs)
+    n_members, n_leadtimes, y, x = pred.shape
+    n_thresholds = len(thr)
+    n_scales = len(scales)
+
+    results = None  # Delay initialization until first score_dict is available
+
+    # Loop over thresholds, scales, members, and lead times
+    for i, thr_val in enumerate(thr):
+        for s, scale in enumerate(scales):
+            for m in range(n_members):
+                for t in range(n_leadtimes):
+                    fcst_bin = (pred[m, t] > thr_val).astype(float)
+                    obs_bin = (obs[t] > thr_val).astype(float)
+                    
+                    # Smooth and re-binarize
+                    S_f, S_o = smooth_fields(fcst_bin, obs_bin, scale)
+                    S_f_bin = S_f > smooth_thr
+                    S_o_bin = S_o > smooth_thr
+                    
+                    # Compute scores using det_cat_fct
+                    score_dict = det_cat_fct(S_f_bin, S_o_bin, thr_val, scores, axis, return_table)
+                    
+                    # Initialize results dictionary on first pass
+                    if results is None:
+                        results = {
+                            key: np.empty((n_thresholds, n_scales, n_members, n_leadtimes))
+                            for key in score_dict.keys()
+                        }
+
+                    # Store each score
+                    for score_name, value in score_dict.items():
+                        results[score_name][i, s, m, t] = value
+
+    return results
+
diff --git a/pysteps/verification/interface.py b/pysteps/verification/interface.py
index 6e325f6f4..d99b5e0c7 100644
--- a/pysteps/verification/interface.py
+++ b/pysteps/verification/interface.py
@@ -87,6 +87,8 @@ def get_method(name, type="deterministic"):
         +------------+--------------------------------------------------------+
         |  FSS       | fractions skill score                                  |
         +------------+--------------------------------------------------------+
+        |  FSS_ens   | fractions skill score for ensemble forecast            |
+        +------------+--------------------------------------------------------+
         |  SAL       | Structure-Amplitude-Location score                     |
         +------------+--------------------------------------------------------+
 
@@ -164,12 +166,12 @@ def get_method(name, type="deterministic"):
     type = type.lower()
 
     if type == "deterministic":
-        from .detcatscores import det_cat_fct
+        from .detcatscores import det_cat_fct, det_cat_fct_spatial
         from .detcontscores import det_cont_fct
-        from .spatialscores import binary_mse, fss
+        from .spatialscores import binary_mse, fss, fss_ens
         from .salscores import sal
 
-        # categorical
+        # categorical or spatial categorical
         if name in [
             "acc",
             "bias",
@@ -184,12 +186,23 @@ def get_method(name, type="deterministic"):
             "pod",
             "sedi",
         ]:
-
             def f(fct, obs, **kwargs):
-                return det_cat_fct(fct, obs, kwargs.pop("thr"), [name])
-
+                if "scales" in kwargs:
+                    return det_cat_fct_spatial(
+                        fct,
+                        obs,
+                        thr=kwargs.pop("thr", 0.0),
+                        scales=kwargs.pop("scales"),
+                        axis=kwargs.pop("axis", None),
+                        scores=[name],
+                        smooth_thr=kwargs.pop("smooth_thr", 0.0),
+                    )[name.upper()]
+                else:
+                    return det_cat_fct(fct, obs, kwargs.pop("thr"), [name])[name.upper()]
+        
             return f
 
+
         # continuous
         elif name in [
             "beta",
@@ -217,6 +230,8 @@ def f(fct, obs, **kwargs):
             return binary_mse
         elif name == "fss":
             return fss
+        elif name == "fss_ens":
+            return fss_ens
         elif name == "sal":
             return sal
 
diff --git a/pysteps/verification/spatialscores.py b/pysteps/verification/spatialscores.py
index 41e02c93e..75a5157d1 100644
--- a/pysteps/verification/spatialscores.py
+++ b/pysteps/verification/spatialscores.py
@@ -677,6 +677,250 @@ def fss_compute(fss):
     return 1.0 - numer / denom
 
 
+def fss_ens(fcst, obs, thr, scales, version='fss'):
+    """
+    Compute Fractions Skill Score (FSS) for ensemble or deterministic forecasts
+    across multiple thresholds and spatial scales.
+
+    Parameters
+    ----------
+    fcst : np.ndarray
+        Forecast field. Shape can be:
+        - (n_members, n_leadtimes, y, x)
+        - (n_leadtimes, y, x)
+        - (n_members, y, x)
+        - (y, x)
+    obs : np.ndarray
+        Observation field. Shape can be:
+        - (n_leadtimes, y, x)
+        - (y, x)
+    thr : float or list of float
+        Threshold(s) used to binarize forecast and observation fields.
+    scales : list or np.ndarray
+        List of spatial scales (window sizes) to apply smoothing.
+    version : str
+        Type of FSS computation:
+        - 'fss'   : member-wise FSS
+        - 'avfss' : average of member-wise FSS
+        - 'pfss'  : FSS of ensemble mean probability
+        - 'agfss' : FSS of aggregated ensemble fields
+
+    Returns
+    -------
+    fss_scores : np.ndarray
+        FSS scores with shape:
+        - (n_thresholds, n_scales, n_members, n_leadtimes) for 'fss'
+        - (n_thresholds, n_scales, 1, n_leadtimes) for ensemble versions
+    """
+    version = version.lower()
+    scales = np.atleast_1d(scales)
+    thr = np.atleast_1d(thr)
+
+    # Normalize input shapes
+    fcst, obs = normalize_forecast_obs(fcst, obs)
+    n_members, n_leadtimes, y, x = fcst.shape
+    n_scales = len(scales)
+    n_thresholds = len(thr)
+
+    # Initialize output array
+    if version == 'fss':
+        fss_scores = np.empty((n_thresholds, n_scales, n_members, n_leadtimes))
+    else:
+        fss_scores = np.empty((n_thresholds, n_scales, 1, n_leadtimes))
+
+    # Loop over thresholds and lead times
+    for i, thr_val in enumerate(thr):
+        for t in range(n_leadtimes):
+            obs_bin = (obs[t] >= thr_val).astype(float)
+
+            if version in ['fss', 'avfss']:
+                for m in range(n_members):
+                    fcst_bin = (fcst[m, t] >= thr_val).astype(float)
+                    for s, scale in enumerate(scales):
+                        num, denom = _compute_fss_components(fcst_bin, obs_bin, scale)
+                        fss_scores[i, s, m, t] = _compute_fss_score(num, denom)
+
+                if version == 'avfss':
+                    for s in range(n_scales):
+                        fss_scores[i, s, 0, t] = fss_scores[i, s, :, t].mean()
+
+            elif version == 'pfss':
+                fcst_bin = (fcst[:, t] >= thr_val).astype(float)
+                fcst_mean = np.mean(fcst_bin, axis=0)
+                for s, scale in enumerate(scales):
+                    num, denom = _compute_fss_components(fcst_mean, obs_bin, scale)
+                    fss_scores[i, s, 0, t] = _compute_fss_score(num, denom)
+
+            elif version == 'agfss':
+                for s, scale in enumerate(scales):
+                    num, denom = 0.0, 0.0
+                    for m in range(n_members):
+                        fcst_bin = (fcst[m, t] >= thr_val).astype(float)
+                        n, d = _compute_fss_components(fcst_bin, obs_bin, scale)
+                        num += n
+                        denom += d
+                    fss_scores[i, s, 0, t] = _compute_fss_score(num, denom)
+
+            else:
+                raise ValueError("Invalid version. Choose from 'fss', 'pfss', 'agfss', 'avfss'.")
+
+    return fss_scores
+
+
+def normalize_forecast_obs(pred, obs):
+    """
+    Normalize forecast and observation arrays to shape:
+    - pred: (n_members, n_leadtimes, y, x)
+    - obs : (n_leadtimes, y, x)
+
+    Parameters
+    ----------
+    pred : np.ndarray
+        Forecast array with shape:
+        - (y, x)
+        - (n_members, y, x)
+        - (n_leadtimes, y, x)
+        - (n_members, n_leadtimes, y, x)
+    obs : np.ndarray
+        Observation array with shape:
+        - (y, x)
+        - (n_leadtimes, y, x)
+
+    Returns
+    -------
+    pred : np.ndarray
+        Normalized forecast array of shape (n_members, n_leadtimes, y, x)
+    obs : np.ndarray
+        Normalized observation array of shape (n_leadtimes, y, x)
+    """
+    pred = np.asarray(pred)
+    obs = np.asarray(obs)
+
+    # Normalize observation shape
+    if obs.ndim == 2:  # (y, x)
+        obs = obs[np.newaxis, ...]
+    elif obs.ndim == 3:
+        pass
+    else:
+        raise ValueError(f"Unsupported observation shape: {obs.shape}")
+
+    n_leadtimes = obs.shape[0]
+
+    # Normalize forecast shape
+    if pred.ndim == 2:  # (y, x)
+        pred = pred[np.newaxis, np.newaxis, ...]
+    elif pred.ndim == 3:
+        if pred.shape[0] == n_leadtimes:  # (n_leadtimes, y, x)
+            pred = pred[np.newaxis, ...]
+        else:  # (n_members, y, x)
+            pred = pred[:, np.newaxis, ...]
+    elif pred.ndim == 4:
+        if pred.shape[1] != n_leadtimes:
+            raise ValueError(f"Forecast lead time dimension {pred.shape[1]} does not match obs {n_leadtimes}")
+    else:
+        raise ValueError(f"Unsupported forecast shape: {pred.shape}")
+
+    return pred, obs
+
+
+def smooth_fields(fcst_bin, obs_bin, scale):
+    """
+    Apply uniform smoothing to forecast and observation binary fields.
+
+    Parameters
+    ----------
+    fcst_bin : np.ndarray
+        Binary forecast field.
+    obs_bin : np.ndarray
+        Binary observation field.
+    scale : int
+        Size of the smoothing window.
+
+    Returns
+    -------
+    S_f : np.ndarray
+        Smoothed forecast field.
+    S_o : np.ndarray
+        Smoothed observation field.
+    """
+    if scale > 1:
+        S_f = uniform_filter(fcst_bin, size=scale, mode="constant", cval=0.0)
+        S_o = uniform_filter(obs_bin, size=scale, mode="constant", cval=0.0)
+    else:
+        S_f = fcst_bin
+        S_o = obs_bin
+    return S_f, S_o
+
+
+def _compute_fss_terms(S_f, S_o):
+    """
+    Compute numerator and denominator for FSS from smoothed fields.
+
+    Parameters
+    ----------
+    S_f : np.ndarray
+        Smoothed forecast field.
+    S_o : np.ndarray
+        Smoothed observation field.
+
+    Returns
+    -------
+    num : float
+        Sum of squared differences.
+    denom : float
+        Sum of squared values.
+    """
+    num = np.sum((S_f - S_o) ** 2)
+    denom = np.sum(S_f ** 2 + S_o ** 2)
+    return num, denom
+
+
+def _compute_fss_components(fcst_bin, obs_bin, scale):
+    """
+    Compute numerator and denominator for the Fractions Skill Score (FSS)
+    at a given spatial scale using smoothed binary fields.
+
+    Parameters
+    ----------
+    fcst_bin : np.ndarray
+        Binary forecast field.
+    obs_bin : np.ndarray
+        Binary observation field.
+    scale : int
+        Size of the smoothing window.
+
+    Returns
+    -------
+    num : float
+        Sum of squared differences.
+    denom : float
+        Sum of squared values.
+    """
+    S_f, S_o = smooth_fields(fcst_bin, obs_bin, scale)
+    return _compute_fss_terms(S_f, S_o)
+
+
+def _compute_fss_score(num, denom, eps=1e-10):
+    """
+    Compute the Fractions Skill Score (FSS) from numerator and denominator.
+
+    Parameters
+    ----------
+    num : float
+        Numerator of the FSS formula (error term).
+    denom : float
+        Denominator of the FSS formula (normalization term).
+    eps : float, optional
+        Small constant to avoid division by zero (default is 1e-10).
+
+    Returns
+    -------
+    fss : float
+        FSS score ranging from 0 (no skill) to 1 (perfect match).
+    """
+    return 1.0 - num / (denom + eps)
+
+
 def _wavelet_decomp(X, w):
     c = pywt.wavedec2(X, w)