Repository: opendatacube/odc-geo
Pull request #143 (2024-04-10 14:59:15)

Fix NaN nodata values in band statistics
When the `save_cog_with_dask` function computed statistics on a float32 band that included NaN values and had a nodata value of NaN, the minimum and maximum values were incorrectly computed as NaN.

When the `save_cog_with_dask` function computed statistics on a float32 band that included NaN values and had a finite nodata value, the NaN values were included in the statistics computation.  Although the correct behavior here is debatable, the best thing to do is to match GDAL and GDAL excludes the NaN values.

This PR updates NaN value handling in statistics computations and adds test coverage to verify how things are handled.

Note that `dask.array`'s `nanmin`, `nanmax`, etc. methods are equivalent to `min`, `max`, etc. when used with integer data.

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/opendatacube/odc-geo/pull/143?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 95.24%. Comparing base [(`6af5d0c`)](https://app.codecov.io/gh/opendatacube/odc-geo/commit/6af5d0cce76d073e83ade6dca51b0377ca7d4e57?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) to head [(`6369525`)](https://app.codecov.io/gh/opendatacube/odc-geo/pull/143?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).
> Report is 5 commits behind head on develop.

> :exclamation: Current head 6369525 differs from pull request most recent head 468da40. Consider uploading reports for the commit 468da40 to get more accurate results


<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop     #143      +/-   ##
===========================================
- Coverage    95.26%   95.24%   -0.02%     
===========================================
  Files           31       31              
  Lines         5323     5326       +3     
===========================================
+ Hits          5071     5073       +2     
- Misses         252      253       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/opendatacube/odc-geo/pull/143?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).


Patch:
diff --git a/odc/geo/cog/_tifffile.py b/odc/geo/cog/_tifffile.py
index 6d95bd9..927c391 100644
--- a/odc/geo/cog/_tifffile.py
+++ b/odc/geo/cog/_tifffile.py
@@ -87,20 +87,8 @@ def _stats_from_layer(
 
     axis = (yaxis, yaxis + 1)
     npix = pix.shape[yaxis] * pix.shape[yaxis + 1]
-    if nodata is not None:
-        dd = da.ma.masked_equal(pix, nodata)
-        return unwrap(
-            {
-                "minimum": dd.min(axis=axis),
-                "maximum": dd.max(axis=axis),
-                "mean": dd.mean(axis=axis),
-                "stddev": dd.std(axis=axis),
-                "valid_percent": da.isfinite(dd).sum(axis=axis) * (100 / npix),
-            },
-            pix.ndim,
-        )
 
-    if pix.dtype.kind == "f":
+    if nodata is None or np.isnan(nodata):
         dd = pix
         return unwrap(
             {
@@ -113,7 +101,8 @@ def _stats_from_layer(
             pix.ndim,
         )
 
-    dd = pix
+    # Exclude both nodata and invalid (e.g. NaN) values from statistics computation
+    dd = da.ma.masked_where((pix == nodata) | ~(np.isfinite(pix)), pix)
     return unwrap(
         {
             "minimum": dd.min(axis=axis),


Test patch:
diff --git a/tests/test_cog.py b/tests/test_cog.py
index ee1561e..378fb34 100644
--- a/tests/test_cog.py
+++ b/tests/test_cog.py
@@ -3,7 +3,9 @@ from io import BytesIO
 from pathlib import Path
 from typing import Optional, Tuple
 
+import numpy as np
 import pytest
+from dask import array as da
 
 from odc.geo.cog import CogMeta, cog_gbox, save_cog_with_dask
 from odc.geo.cog._shared import compute_cog_spec, num_overviews
@@ -14,6 +16,7 @@ from odc.geo.cog._tifffile import (
     _gdal_sample_descriptions,
     _make_empty_cog,
     _norm_compression_tifffile,
+    _stats_from_layer,
     geotiff_metadata,
 )
 from odc.geo.geobox import GeoBox
@@ -450,3 +453,57 @@ def test_cog_with_dask_smoke_test(gbox: GeoBox, tmp_path: Path, dtype):
     fname = str(tmp_path / "cog-bandnames-incorrect.tif")
     with pytest.raises(ValueError):
         save_cog_with_dask(img, fname, compression="deflate", level=2)
+
+
+@pytest.mark.parametrize(
+    ("array", "nodata", "minimum", "maximum", "mean", "stddev", "valid_percent"),
+    [
+        pytest.param([[1, 1], [1, 1]], None, 1, 1, 1, 0, 100, id="basic int"),
+        pytest.param(
+            [[1.0, 1.0], [1.0, 1.0]], None, 1.0, 1.0, 1.0, 0.0, 100, id="basic float"
+        ),
+        pytest.param([[1, 0], [0, 1]], 0, 1, 1, 1, 0, 50, id="int with numeric nodata"),
+        pytest.param(
+            [[1.0, 0.0], [0.0, 1.0]], 0, 1, 1, 1, 0, 50, id="float with numeric nodata"
+        ),
+        pytest.param(
+            [[1.0, np.nan], [np.nan, 1.0]],
+            None,
+            1,
+            1,
+            1,
+            0,
+            50,
+            id="float with nan, None nodata",
+        ),
+        pytest.param(
+            [[1.0, np.nan], [np.nan, 1.0]],
+            np.nan,
+            1,
+            1,
+            1,
+            0,
+            50,
+            id="float with nan, nan nodata",
+        ),
+        pytest.param(
+            [[1.0, np.nan], [0, 1.0]],
+            0,
+            1,
+            1,
+            1,
+            0,
+            50,
+            id="float with nan, numeric nodata",
+        ),
+    ],
+)
+def test_stats_from_layer(array, nodata, minimum, maximum, mean, stddev, valid_percent):
+    x = da.from_array(array)
+    stats = _stats_from_layer(x, nodata).compute()[0]
+
+    assert stats["minimum"] == minimum
+    assert stats["maximum"] == minimum
+    assert stats["mean"] == mean
+    assert stats["stddev"] == stddev
+    assert stats["valid_percent"] == valid_percent