Add tests and fix block-level M2 for variance stability (#1090)

- Block-level sum_squares now computes M2 (sum of squared deviations
  from block mean) instead of raw sum(x²), avoiding float64 precision
  loss for large values.
- Updated test_stats_all_nan_zone and test_stats_nodata_wipes_zone to
  expect NaN from dask (no longer 0).
- Added test_stats_variance_numerical_stability_1090: values near 1e8
  with spread of 1, verifying dask matches numpy to 1e-6.
- Added test_stats_nodata_none_no_warning_1090: confirms no
  FutureWarning when nodata_values=None.

Author: brendancol

xrspatial/tests/test_zonal.py

@@ -680,18 +680,8 @@ def test_stats_all_nan_zone(backend):
             'sum':   [12.0],
             'count': [2],
         }
-    elif 'dask' in backend:
-        # dask uses nansum reduction, so count/sum of all-NaN become 0
-        expected = {
-            'zone':  [1, 2],
-            'mean':  [np.nan, 6.0],
-            'max':   [np.nan, 7.0],
-            'min':   [np.nan, 5.0],
-            'sum':   [0.0, 12.0],
-            'count': [0, 2],
-        }
     else:
-        # numpy keeps empty zone with NaN for every stat
+        # numpy and dask both return NaN for all-NaN zones
         expected = {
             'zone':  [1, 2],
             'mean':  [np.nan, 6.0],
@@ -798,16 +788,8 @@ def test_stats_nodata_wipes_zone(backend):
             'sum':   [10.0],
             'count': [2],
         }
-    elif 'dask' in backend:
-        expected = {
-            'zone':  [1, 2],
-            'mean':  [np.nan, 5.0],
-            'max':   [np.nan, 7.0],
-            'min':   [np.nan, 3.0],
-            'sum':   [0.0, 10.0],
-            'count': [0, 2],
-        }
     else:
+        # numpy and dask both return NaN for zones with no valid values
         expected = {
             'zone':  [1, 2],
             'mean':  [np.nan, 5.0],
@@ -868,6 +850,71 @@ def test_zonal_stats_inputs_unmodified(backend, data_zones, data_values_2d, resu
     assert_input_data_unmodified(data_values_2d, copied_data_values_2d)
 
 
+@pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
+@pytest.mark.filterwarnings("ignore:invalid value encountered in divide:RuntimeWarning")
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+def test_stats_variance_numerical_stability_1090(backend):
+    """Dask std/var should match numpy for data with large mean, small spread.
+
+    Regression test for #1090: the naive one-pass formula
+    ``(Σx² − (Σx)²/n) / n`` loses precision through catastrophic
+    cancellation.  The fix uses Chan-Golub-LeVeque parallel merge.
+    """
+    if 'dask' in backend and not dask_array_available():
+        pytest.skip("Requires Dask")
+
+    # Values near 1e8 with a spread of 1: the naive formula would lose
+    # most of the significant digits in float64.
+    zones_data = np.array([[1, 1, 1, 1, 1, 1]])
+    values_data = np.array([[1e8, 1e8 + 1, 1e8 + 2,
+                             1e8 + 3, 1e8 + 4, 1e8 + 5]], dtype=np.float64)
+
+    zones = create_test_raster(zones_data, backend, chunks=(1, 3))
+    values = create_test_raster(values_data, backend, chunks=(1, 3))
+
+    df_result = stats(zones=zones, values=values,
+                      stats_funcs=['mean', 'std', 'var'])
+
+    if hasattr(df_result, 'compute'):
+        df_result = df_result.compute()
+
+    # Reference: population variance of [0,1,2,3,4,5] = 35/12 ≈ 2.9167
+    expected_var = np.var(np.arange(6, dtype=np.float64))
+    expected_std = np.std(np.arange(6, dtype=np.float64))
+
+    actual_var = float(df_result['var'].iloc[0])
+    actual_std = float(df_result['std'].iloc[0])
+
+    assert abs(actual_var - expected_var) < 1e-6, (
+        f"var={actual_var}, expected={expected_var}"
+    )
+    assert abs(actual_std - expected_std) < 1e-6, (
+        f"std={actual_std}, expected={expected_std}"
+    )
+
+
+def test_stats_nodata_none_no_warning_1090():
+    """Passing nodata_values=None (the default) should not trigger warnings.
+
+    Regression test for #1090: ``zone_values != None`` triggered a numpy
+    FutureWarning.
+    """
+    import warnings
+
+    zones_data = np.array([[1, 1], [2, 2]], dtype=float)
+    values_data = np.array([[1.0, 2.0], [3.0, 4.0]])
+    zones = xr.DataArray(zones_data)
+    values = xr.DataArray(values_data)
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        df = stats(zones=zones, values=values, nodata_values=None)
+
+    assert len(df) == 2
+    assert float(df['mean'].iloc[0]) == 1.5
+    assert float(df['mean'].iloc[1]) == 3.5
+
+
 @pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
 @pytest.mark.filterwarnings("ignore:invalid value encountered in divide:RuntimeWarning")
 @pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])

xrspatial/zonal.py

@@ -198,7 +198,7 @@ def _stats_majority(data):
     min=lambda z: z.min(),
     sum=lambda z: z.sum(),
     count=lambda z: _stats_count(z),
-    sum_squares=lambda z: (z**2).sum()
+    sum_squares=lambda z: ((z - z.mean()) ** 2).sum()  # block-level M2
 )
 
 
@@ -227,15 +227,13 @@ def _dask_mean(sums, counts):  # noqa
     return sums / counts
 
 
-def _parallel_variance(block_counts, block_sums, block_sum_squares):
+def _parallel_variance(block_counts, block_sums, block_m2s):
     """Population variance via Chan-Golub-LeVeque parallel merge.
 
-    Each input is (n_blocks, n_zones).  Returns (n_zones,) variance,
+    Each input is (n_blocks, n_zones).  ``block_m2s`` contains
+    per-block M2 values (sum of squared deviations from the block mean),
+    NOT raw sum-of-squares.  Returns (n_zones,) population variance,
     with NaN for zones that have no valid values in any block.
-
-    This avoids the naive ``(Σx² − (Σx)²/n) / n`` formula whose
-    subtraction can lose most significant digits when the mean is
-    large relative to the standard deviation.
     """
     n_blocks = block_counts.shape[0]
     n_zones = block_counts.shape[1]
@@ -247,14 +245,13 @@ def _parallel_variance(block_counts, block_sums, block_sum_squares):
     for i in range(n_blocks):
         nc = np.asarray(block_counts[i], dtype=np.float64)
         sc = np.asarray(block_sums[i], dtype=np.float64)
-        sqc = np.asarray(block_sum_squares[i], dtype=np.float64)
+        m2_b = np.asarray(block_m2s[i], dtype=np.float64)
 
         has_data = np.isfinite(nc) & (nc > 0)
         nc_safe = np.where(has_data, nc, 1.0)  # avoid /0
 
         with np.errstate(invalid='ignore', divide='ignore'):
             mean_b = sc / nc_safe
-            m2_b = sqc - sc ** 2 / nc_safe  # block-internal M2
 
         nc = np.where(has_data, nc, 0.0)
         n_ab = n_acc + nc