Repository: opendatacube/datacube-core

Issue:
Numpy>2.0 causes loss of precision during reproject_and_fuse for inputs with very small resolutions. 
### Background 

We have some very high resolution rasters (4.2141938e-08 resolution) in EPSG:4326 that were able to be reprojected/resampled/etc. prior to numpy 2.0.


### Expected behaviour

A correct reprojection/resampling operation for input sources with small resolution values.


### Actual behaviour

numpy 2.0+ started preserving dtypes when doing scalar operations [ref](https://numpy.org/devdocs/numpy_2_0_migration_guide.html#changes-to-numpy-data-type-promotion). The [`gbox_boundary`](https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/geometry/tools.py#L58) func produces float32 coordinates that were previously implicitly converted to float64 when [computing the roi](https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/geometry/tools.py#L617). With numpy>2.0 this computation is done on the source float32 array and the pixel space coordinates being produced are no longer correct.  


### Steps to reproduce the behaviour

A minimal reproducible example is below. Numpy 2.1.3 was used for this example.
```
import numpy as np

from datacube.utils.geometry._base import GeoBox, CRS, Affine
from datacube.utils.geometry.tools import (
    compute_reproject_roi,
    native_pix_transform,
    gbox_boundary,
    roi_from_points,
    roi_is_empty,
    roi_boundary,
)


src_gbox = GeoBox(
    4096,
    4096,
    Affine(
        4.214193778650287e-08,
        0.0,
        -52.78332070428918,
        0.0,
        -4.2139759941471435e-08,
        -24.400836166256216,
    ),
    CRS(
        'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
    ),
)
dst_gbox = GeoBox(
    2390,
    2390,
    Affine(0.001, 0.0, -5875812.3, 0.0, -0.001, -2802348.662),
    CRS("EPSG:3857"),
)

# defining some constants from the func defaults
pts_per_side = 5

tr = native_pix_transform(src_gbox, dst_gbox)

# gbox_boundary(dst_gbox, pts_per_side) doesn't change with the np version
tr.back(gbox_boundary(dst_gbox, pts_per_side)) # this version is incorrect and results in gaps in the output dataset
tr.back(gbox_boundary(dst_gbox, pts_per_side).astype("float64")) # this version produces correct results
```

With numpy 2.1.3, the output of tr.back is:
```
[(-24.901943922042847, 3956.8316799402237),
 (81.68016076087952, 3956.8316799402237),
 (188.2622652053833, 3956.8316799402237),
 (401.426474571228, 3956.8316799402237),
 (508.0085790157318, 3956.8316799402237),
 (508.0085790157318, 4053.898505091667),
 (508.0085790157318, 4199.498736977577),
 (508.0085790157318, 4296.565554261208),
 (508.0085790157318, 4442.165774464607),
 (401.426474571228, 4442.165774464607),
 (188.2622652053833, 4442.165774464607),
 (81.68016076087952, 4442.165774464607),
 (-24.901943922042847, 4442.165774464607),
 (-24.901943922042847, 4296.565554261208),
 (-24.901943922042847, 4199.498736977577),
 (-24.901943922042847, 4053.898505091667)]
```
vs either numpy<2.0 or with the float64 cast:
```
[(17.730897903442383, 3939.7479182481766),
 (145.0965130329132, 3939.7479182481766),
 (272.46212792396545, 3939.7479182481766),
 (399.8277428150177, 3939.7479182481766),
 (527.19335770607, 3939.7479182481766),
 (527.19335770607, 4055.7427748441696),
 (527.19335770607, 4171.737626671791),
 (527.19335770607, 4287.732474088669),
 (527.19335770607, 4403.727317094803),
 (399.8277428150177, 4403.727317094803),
 (272.46212792396545, 4403.727317094803),
 (145.0965130329132, 4403.727317094803),
 (17.730897903442383, 4403.727317094803),
 (17.730897903442383, 4287.732474088669),
 (17.730897903442383, 4171.737626671791),
 (17.730897903442383, 4055.7427748441696)]
```

### Environment information

datacube==1.8.19
pyproj==3.7.0
numpy versions mentioned above are 1.26.4 and 2.1.3


### Potential solutions

Updating `gbox_boundary`/`roi_boundary` to create float64 arrays would solve this - since it's only creating sets of 16 points per call I don't think this would cause significant performance impacts.  

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1673?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 [(`cd5f4e9`)](https://app.codecov.io/gh/opendatacube/datacube-core/commit/cd5f4e99d181718050b8563364f02ee2424e4a72?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube) to head [(`8438319`)](https://app.codecov.io/gh/opendatacube/datacube-core/commit/843831996f19a2cbaea907a36e176f211a11e98a?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=opendatacube).

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


```diff
@@           Coverage Diff            @@
##           develop    #1673   +/-   ##
========================================
  Coverage    91.78%   91.78%           
========================================
  Files          132      132           
  Lines        14646    14646           
========================================
  Hits         13443    13443           
  Misses        1203     1203           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/opendatacube/datacube-core/pull/1673?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).

SpacemanPaul: No probs with this, but needs whats_new entry.

Patch:
diff --git a/datacube/utils/geometry/tools.py b/datacube/utils/geometry/tools.py
index c4c7018e..8b352245 100644
--- a/datacube/utils/geometry/tools.py
+++ b/datacube/utils/geometry/tools.py
@@ -60,8 +60,8 @@ def gbox_boundary(gbox, pts_per_side=16):
 
     """
     H, W = gbox.shape[:2]
-    xx = np.linspace(0, W, pts_per_side, dtype='float32')
-    yy = np.linspace(0, H, pts_per_side, dtype='float32')
+    xx = np.linspace(0, W, pts_per_side, dtype='float64')
+    yy = np.linspace(0, H, pts_per_side, dtype='float64')
 
     return polygon_path(xx, yy).T[:-1]
 
@@ -72,11 +72,11 @@ def roi_boundary(roi, pts_per_side=2):
 
     roi needs to be in the normalised form, i.e. no open-ended start/stop, see roi_normalise
 
-    :returns: Nx2 float32 array of X,Y points on the perimeter of the envelope defined by `roi`
+    :returns: Nx2 float64 array of X,Y points on the perimeter of the envelope defined by `roi`
     """
     yy, xx = roi
-    xx = np.linspace(xx.start, xx.stop, pts_per_side, dtype='float32')
-    yy = np.linspace(yy.start, yy.stop, pts_per_side, dtype='float32')
+    xx = np.linspace(xx.start, xx.stop, pts_per_side, dtype='float64')
+    yy = np.linspace(yy.start, yy.stop, pts_per_side, dtype='float64')
 
     return polygon_path(xx, yy).T[:-1]
 
diff --git a/docs/about/whats_new.rst b/docs/about/whats_new.rst
index d10b2bff..7992dedd 100644
--- a/docs/about/whats_new.rst
+++ b/docs/about/whats_new.rst
@@ -11,6 +11,7 @@ v1.8.next
 - Update docs (:pull:`1631`,:pull:`1651`)
 - Update docker image and CI (multiple PRs)
 - Pin out bad versions of dask (:pull:`1663`)
+- Ensure correct boundary scalar precision with numpy 2 (:pull:`1673`)
 
 v1.8.19 (2nd July 2024)
 =======================
diff --git a/wordlist.txt b/wordlist.txt
index be3c5e7e..50d6cc46 100644
--- a/wordlist.txt
+++ b/wordlist.txt
@@ -377,6 +377,7 @@ PRIMEM
 prog
 Proj
 provence
+PRs
 psql
 pts
 pwd


Test patch:
diff --git a/tests/test_geometry.py b/tests/test_geometry.py
index 5c5721b9..2382c987 100644
--- a/tests/test_geometry.py
+++ b/tests/test_geometry.py
@@ -50,6 +50,7 @@ from datacube.utils.geometry._base import (
     _norm_crs_or_error,
     _make_crs_key,
 )
+from datacube.utils.geometry.tools import gbox_boundary
 from datacube.testutils.geom import (
     epsg4326,
     epsg3577,
@@ -1588,3 +1589,50 @@ def test_lonalt_bounds_more_than_180():
 
     assert geometry.lonlat_bounds(poly, "quick") == approx((-150, -30, 150, 30))
     assert geometry.lonlat_bounds(poly, "safe") == approx((-150, -30, 150, 30))
+
+
+def test_small_resolution_precision():
+    src_gbox = GeoBox(
+        4096,
+        4096,
+        Affine(
+            4.214193778650287e-08,
+            0.0,
+            -52.78332070428918,
+            0.0,
+            -4.2139759941471435e-08,
+            -24.400836166256216,
+        ),
+        CRS(
+            """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,
+            AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
+            UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],
+            AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"""
+        ),
+    )
+    dst_gbox = GeoBox(
+        2390,
+        2390,
+        Affine(0.001, 0.0, -5875812.3, 0.0, -0.001, -2802348.662),
+        CRS("EPSG:3857"),
+    )
+
+    tr = native_pix_transform(src_gbox, dst_gbox)
+    assert tr.back(gbox_boundary(dst_gbox, 5)) == [
+        (17.730897903442383, 3939.7479182481766),
+        (145.0965130329132, 3939.7479182481766),
+        (272.46212792396545, 3939.7479182481766),
+        (399.8277428150177, 3939.7479182481766),
+        (527.19335770607, 3939.7479182481766),
+        (527.19335770607, 4055.7427748441696),
+        (527.19335770607, 4171.737626671791),
+        (527.19335770607, 4287.732474088669),
+        (527.19335770607, 4403.727317094803),
+        (399.8277428150177, 4403.727317094803),
+        (272.46212792396545, 4403.727317094803),
+        (145.0965130329132, 4403.727317094803),
+        (17.730897903442383, 4403.727317094803),
+        (17.730897903442383, 4287.732474088669),
+        (17.730897903442383, 4171.737626671791),
+        (17.730897903442383, 4055.7427748441696)
+    ]