Repository: opendatacube/datacube-core

Issue:
Important fix proposal: Implement setting "nan" nodata value for floating point rasters with no nodata attribute 
Hi all, as discussed internally within DEA, we have recently encountered an issue where floating point rasters are not displayed correctly in ESRI software due to the lack of an explicitly set nodata attribute of "nan". This differs from open software such as QGIS and GDAL, which correctly interpret nodata by implicitly treating a missing nodata attribute as equivalent to "nan" (**edit:** this turns out to not be the case: [GDAL overview generation is also affected by this issue](https://github.com/opendatacube/datacube-core/issues/1595#issuecomment-2167274760)).

While this issue has been resolved by patching previously generated data on prod S3, a more sustainable long-term fix would be to set nodata attributes in our GeoTIFF writing code, to ensure that any floating point raster without a custom nodata attribute is assigned a nodata attribute of "nan".

The `datacube.utils.cog.write_cog` and `datacube.utils.cog.to_cog` functions are important places to fix this, as they are used in the generation of almost all ODC products. A possible fix would likely involve a minor change similar to psuedocode below [somewhere here](https://github.com/opendatacube/datacube-core/blob/8b46c2284e13c74739ca16faba67669af0b9cf74/datacube/utils/cog.py#L210-L362):

```
if (raster is floating point dtype) AND (no nodata attribute set):
    raster.nodata = "nan"
```

The alternative to fixing it in GeoTIFF writing code would be to make sure all downstream product worklows manually set nodata values. However, in my opinion this would be difficult to document and would be highly likely to be missed in the future, leading to more expensive remediation. Fixing it at the source seems like a better option, especially as extensive testing has revealed no downsteam impacts of the change while successfully solve the issue for ESRI users (a large proportion of our user base). 

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1602?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 91.78%. Comparing base [(`8b46c22`)](https://app.codecov.io/gh/opendatacube/datacube-core/commit/8b46c2284e13c74739ca16faba67669af0b9cf74?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) to head [(`f01146e`)](https://app.codecov.io/gh/opendatacube/datacube-core/commit/f01146e026622dfa316d0601b66768e5d3d09b63?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).

> :exclamation: **Current head f01146e differs from pull request most recent head 783124f**
> 
> Please [upload](https://docs.codecov.com/docs/codecov-uploader) reports for the commit 783124f to get more accurate results.


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


```diff
@@           Coverage Diff            @@
##           develop    #1602   +/-   ##
========================================
  Coverage    91.78%   91.78%           
========================================
  Files          132      132           
  Lines        14637    14640    +3     
========================================
+ Hits         13434    13437    +3     
  Misses        1203     1203           
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1602?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/.github/workflows/main.yml b/.github/workflows/main.yml
index 69235df1..37a245a0 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -76,7 +76,7 @@ jobs:
         password: ${{ secrets.DOCKERHUBPASSWD }}
 
     - name: Build Docker
-      uses: docker/build-push-action@v5
+      uses: docker/build-push-action@v6
       with:
         file: docker/Dockerfile
         context: .
@@ -108,7 +108,7 @@ jobs:
         github.event_name == 'push'
         && github.ref == 'refs/heads/develop'
         && steps.changes.outputs.docker == 'true'
-      uses: docker/build-push-action@v5
+      uses: docker/build-push-action@v6
       with:
         file: docker/Dockerfile
         context: .
diff --git a/datacube/utils/cog.py b/datacube/utils/cog.py
index c3611eff..5e1b6615 100644
--- a/datacube/utils/cog.py
+++ b/datacube/utils/cog.py
@@ -133,6 +133,10 @@ def _write_cog(
         compress="DEFLATE",
     )
 
+    # If nodata is not set, but the array is of floating point type, force nodata=nan
+    if nodata is None and np.issubdtype(pix.dtype, np.floating):
+        nodata = np.nan
+
     if nodata is not None:
         rio_opts.update(nodata=nodata)
 
diff --git a/datacube/utils/math.py b/datacube/utils/math.py
index 8181bd2c..1d62f229 100644
--- a/datacube/utils/math.py
+++ b/datacube/utils/math.py
@@ -178,6 +178,7 @@ def num2numpy(x, dtype, ignore_range=None):
     :param x int|float: Numerical value to convert to numpy.type
     :param dtype str|numpy.dtype|numpy.type: Destination dtype
     :param ignore_range: If set to True skip range check and cast anyway (for example: -1 -> 255)
+                         (Not supported in numpy 2.0+)
 
     :returns: None if x is None
     :returns: None if x is outside the valid range of dtype and ignore_range is not set
diff --git a/docs/about/whats_new.rst b/docs/about/whats_new.rst
index 9edafd9f..a6366449 100644
--- a/docs/about/whats_new.rst
+++ b/docs/about/whats_new.rst
@@ -8,6 +8,7 @@ What's New
 v1.8.next
 =========
 
+- Always write floating point bands to cogs with nodata=nan for ESRI and GDAL compatibility (:pull:`1602`)
 - Add deprecation warning for config environment names that will not be supported in 1.9 (:pull:`1592`)
 - Update docker image to GDAL 3.9/Python 3.12/Ubuntu 24.04 (:pull:`1587`)
 - Update readthedocs stylesheet for dark theme (:pull:`1579`)
diff --git a/setup.py b/setup.py
index 8464f799..43ee11fb 100755
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@ tests_require = [
     'pytest-cov',
     'pytest-timeout',
     'pytest-httpserver',
-    'moto',
+    'moto<5.0',  # 5.0 will require changes to some tests.
 ]
 doc_require = [
     'Sphinx',
diff --git a/wordlist.txt b/wordlist.txt
index aeeea4c7..afa81ed6 100644
--- a/wordlist.txt
+++ b/wordlist.txt
@@ -150,6 +150,7 @@ EP
 EPSG
 epsg
 ESPA
+ESRI
 evironments
 f'file
 fd


Test patch:
diff --git a/datacube/testutils/__init__.py b/datacube/testutils/__init__.py
index 29637521..91755a88 100644
--- a/datacube/testutils/__init__.py
+++ b/datacube/testutils/__init__.py
@@ -415,11 +415,12 @@ def gen_tiff_dataset(bands,
                            **kwargs)
 
         gbox = meta.gbox
-
-        mm.append(dict(name=name,
-                       path=fname,
-                       layer=1,
-                       dtype=meta.dtype))
+        bdict = dict(name=name,
+                     path=fname,
+                     layer=1,
+                     nodata=band.nodata,
+                     dtype=meta.dtype)
+        mm.append(bdict)
 
     uri = Path(base_folder_of_record/'metadata.yaml').absolute().as_uri()
     ds = mk_sample_dataset(mm,
diff --git a/tests/test_3d.py b/tests/test_3d.py
index 072aa8d9..90b9e979 100644
--- a/tests/test_3d.py
+++ b/tests/test_3d.py
@@ -68,8 +68,8 @@ def test_extra_dimensions(eo3_metadata, cover_z_dataset_type):
     # Check chunk size
     assert dt.extra_dimensions.chunk_size() == (("z",), (30,))
 
-    # String representation
-    readable = (
+    # String representation (numpy 1.x and numpy 2.x formats - numpy 2.x reports total size in bytes)
+    readable_1 = (
         "ExtraDimensions(extra_dim={'z': {'name': 'z', 'values': [5, 10, 15, "
         "20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, "
         "105, 110, 115, 120, 125, 130, 135, 140, 145, 150], 'dtype': "
@@ -80,8 +80,19 @@ def test_extra_dimensions(eo3_metadata, cover_z_dataset_type):
         "140., 145., 150.])\nCoordinates:\n  * z        (z) int64 5 10 15 20 "
         "25 30 35 40 ... 120 125 130 135 140 145 150} )"
     )
-    assert str(dt.extra_dimensions) == readable
-    assert f"{dt.extra_dimensions!r}" == readable
+    readable_2 = (
+        "ExtraDimensions(extra_dim={'z': {'name': 'z', 'values': [5, 10, 15, "
+        "20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, "
+        "105, 110, 115, 120, 125, 130, 135, 140, 145, 150], 'dtype': "
+        "'float64'}}, dim_slice={'z': (0, 30)} coords={'z': <xarray.DataArray "
+        "'z' (z: 30)> Size: 240B\narray([  5.,  10.,  15.,  20.,  25.,  30.,  35.,  "
+        "40.,  45.,  50.,  55.,\n        60.,  65.,  70.,  75.,  80.,  85.,  "
+        "90.,  95., 100., 105., 110.,\n       115., 120., 125., 130., 135., "
+        "140., 145., 150.])\nCoordinates:\n  * z        (z) int64 240B 5 10 15 20 "
+        "25 30 35 ... 120 125 130 135 140 145 150} )"
+    )
+    assert str(dt.extra_dimensions) in (readable_1, readable_2)
+    assert f"{dt.extra_dimensions!r}" in (readable_1, readable_2)
 
 
 def test_extra_dimensions_exceptions(eo3_metadata, cover_z_dataset_type):
diff --git a/tests/test_utils_cog.py b/tests/test_utils_cog.py
index d73e33fb..5b2ab81f 100644
--- a/tests/test_utils_cog.py
+++ b/tests/test_utils_cog.py
@@ -3,6 +3,7 @@
 # Copyright (c) 2015-2024 ODC Contributors
 # SPDX-License-Identifier: Apache-2.0
 import pytest
+import math
 from pathlib import Path
 import numpy as np
 import xarray as xr
@@ -19,8 +20,8 @@ from datacube.testutils.io import native_load, rio_slurp_xarray, rio_slurp
 from datacube.utils.cog import write_cog, to_cog, _write_cog
 
 
-def gen_test_data(prefix, dask=False, shape=None):
-    w, h, dtype, nodata, ndw = 96, 64, "int16", -999, 7
+def gen_test_data(prefix, dask=False, shape=None, dtype="int16", nodata=-999):
+    w, h, ndw = 96, 64, 7
     if shape is not None:
         h, w = shape
 
@@ -109,6 +110,20 @@ def test_cog_file(tmpdir, opts):
     with pytest.warns(UserWarning):
         write_cog(xx, pp / "cog_badblocksize.tif", blocksize=50)
 
+    # check writing floating point COG with no explicit nodata
+    zz, ds = gen_test_data(pp, dtype="float32", nodata=None)
+    # write to file
+    ff = write_cog(
+        zz,
+        pp / "cog_float.tif",
+        **opts
+    )
+    assert isinstance(ff, Path)
+    assert ff == pp / "cog_float.tif"
+    assert ff.exists()
+    aa = rio_slurp_xarray(pp / "cog_float.tif")
+    assert aa.attrs["nodata"] == "nan" or math.isnan(aa.attrs["nodata"])
+
 
 def test_cog_file_dask(tmpdir):
     pp = Path(str(tmpdir))
diff --git a/tests/test_utils_other.py b/tests/test_utils_other.py
index 60ca28fa..462a7553 100644
--- a/tests/test_utils_other.py
+++ b/tests/test_utils_other.py
@@ -588,7 +588,12 @@ def test_num2numpy():
     assert num2numpy(256, 'uint8') is None
     assert num2numpy(-1, 'uint16') is None
     assert num2numpy(-1, 'uint32') is None
-    assert num2numpy(-1, 'uint8', ignore_range=True) == np.uint8(255)
+    try:
+        # Numpy 1.x supports wrapping of unsisinged types
+        assert num2numpy(-1, 'uint8', ignore_range=True) == np.uint8(255)
+    except OverflowError:
+        # Numpy 2.0 will throw Overflow error rather than wrapping
+        pass
 
     assert num2numpy(0, 'uint8') == 0
     assert num2numpy(255, 'uint8') == 255