diff --git a/xrspatial/tests/test_classify.py b/xrspatial/tests/test_classify.py
index d563d9685..651af4be7 100644
--- a/xrspatial/tests/test_classify.py
+++ b/xrspatial/tests/test_classify.py
@@ -1252,3 +1252,102 @@ def test_generate_sample_indices_large_is_deterministic():
     a = _generate_sample_indices(20_000_000, 20_000)
     b = _generate_sample_indices(20_000_000, 20_000)
     np.testing.assert_array_equal(a, b)
+
+
+# ===================================================================
+# Degenerate all-NaN input for the remaining classifiers
+# ===================================================================
+# equal_interval and natural_breaks already cover all-NaN. The other
+# classifiers were not exercised on an all-non-finite raster, where the
+# finite mask removes every element. std_mean and maximum_breaks degrade
+# to an all-NaN result; on the eager (numpy/cupy) backends head_tail_breaks,
+# percentiles, and box_plot currently raise an opaque reduction error
+# (issue #3510), so their tests are xfail until that is fixed. Flip them to
+# plain assertions when #3510 lands. strict=False so a concurrent fix does
+# not break main via XPASS. See the dask section below for the per-backend
+# split (head_tail_breaks already degrades cleanly on dask).
+
+def test_std_mean_all_nan():
+    import warnings
+    agg = xr.DataArray(np.full((4, 5), np.nan))
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore', RuntimeWarning)
+        result = std_mean(agg)
+    assert np.all(np.isnan(result.data))
+
+
+def test_maximum_breaks_all_nan():
+    agg = xr.DataArray(np.full((4, 5), np.nan))
+    result = maximum_breaks(agg, k=5)
+    assert np.all(np.isnan(result.data))
+
+
+@pytest.mark.xfail(reason="all-NaN input crashes; see issue #3510", strict=False)
+def test_head_tail_breaks_all_nan():
+    agg = xr.DataArray(np.full((4, 5), np.nan))
+    result = head_tail_breaks(agg)
+    assert np.all(np.isnan(result.data))
+
+
+@pytest.mark.xfail(reason="all-NaN input crashes; see issue #3510", strict=False)
+def test_percentiles_all_nan():
+    agg = xr.DataArray(np.full((4, 5), np.nan))
+    result = percentiles(agg)
+    assert np.all(np.isnan(result.data))
+
+
+@pytest.mark.xfail(reason="all-NaN input crashes; see issue #3510", strict=False)
+def test_box_plot_all_nan():
+    agg = xr.DataArray(np.full((4, 5), np.nan))
+    result = box_plot(agg)
+    assert np.all(np.isnan(result.data))
+
+
+# All-NaN on the dask backend. The dask paths are separate implementations,
+# and they do not match the eager paths on this degenerate input:
+# std_mean, maximum_breaks, and head_tail_breaks all return all-NaN on dask
+# (head_tail_breaks's dask path has a total_count == 0 guard the eager path
+# lacks, so it does not hit the #3510 crash), while percentiles and box_plot
+# crash on dask too. Pin all of it so the per-backend behaviour is explicit
+# and the #3510 fix can target only the eager paths.
+
+@dask_array_available
+def test_std_mean_all_nan_dask():
+    import warnings
+    agg = xr.DataArray(da.full((4, 5), np.nan, chunks=(2, 5)))
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore', RuntimeWarning)
+        result = std_mean(agg)
+    assert np.all(np.isnan(result.data.compute()))
+
+
+@dask_array_available
+def test_maximum_breaks_all_nan_dask():
+    agg = xr.DataArray(da.full((4, 5), np.nan, chunks=(2, 5)))
+    result = maximum_breaks(agg, k=5)
+    assert np.all(np.isnan(result.data.compute()))
+
+
+@dask_array_available
+def test_head_tail_breaks_all_nan_dask():
+    # Diverges from the eager path: dask returns all-NaN cleanly (#3510 is
+    # eager-only). Plain passing test, not xfail.
+    agg = xr.DataArray(da.full((4, 5), np.nan, chunks=(2, 5)))
+    result = head_tail_breaks(agg)
+    assert np.all(np.isnan(result.data.compute()))
+
+
+@dask_array_available
+@pytest.mark.xfail(reason="all-NaN input crashes; see issue #3510", strict=False)
+def test_percentiles_all_nan_dask():
+    agg = xr.DataArray(da.full((4, 5), np.nan, chunks=(2, 5)))
+    result = percentiles(agg)
+    assert np.all(np.isnan(result.data.compute()))
+
+
+@dask_array_available
+@pytest.mark.xfail(reason="all-NaN input crashes; see issue #3510", strict=False)
+def test_box_plot_all_nan_dask():
+    agg = xr.DataArray(da.full((4, 5), np.nan, chunks=(2, 5)))
+    result = box_plot(agg)
+    assert np.all(np.isnan(result.data.compute()))