Round gridspec tilesize to fix floating point precision issue and related fixes

Author: alexgleith

.github/workflows/main.yml

@@ -78,7 +78,7 @@ jobs:
         run: |
           conda info
           conda list
-          mamba -V
+          # mamba -V
           conda config --show-sources
           conda config --show
           printenv | sort

odc/geo/_dask.py

@@ -1,5 +1,5 @@
 from functools import partial
-from typing import Any, Dict, Optional, Sequence, Tuple, Union
+from typing import Any, Dict, Optional, Tuple, Union
 from uuid import uuid4
 
 import dask.array as da
@@ -15,7 +15,7 @@
 
 
 def _do_chunked_reproject(
-    d2s: Dict[Tuple[int, int], Sequence[Tuple[int, int]]],
+    d2s: Dict[Tuple[int, int], list[Tuple[int, int]]],
     src_gbt: GeoboxTiles,
     dst_gbt: GeoboxTiles,
     dst_idx: Tuple[int, int],

odc/geo/gridspec.py

@@ -66,9 +66,12 @@ def __init__(
         self.crs = norm_crs_or_error(crs)
         self._shape = tile_shape
         self.resolution = resolution
+
+        # This is an arbitrary rounding of 12 decimal places
+        # I don't know how to make it more robust
         self.tile_size = xy_(
-            tile_shape.x * abs(resolution.x),
-            tile_shape.y * abs(resolution.y),
+            round(tile_shape.x * abs(resolution.x), 12),
+            round(tile_shape.y * abs(resolution.y), 12),
         )
         self.origin = origin
 

tests/test_gridspec.py

@@ -9,7 +9,7 @@
 import pytest
 from pytest import approx
 
-from odc.geo import CRS, res_, resyx_, xy_, yx_
+from odc.geo import CRS, res_, resyx_, xy_, yx_, XY
 from odc.geo.geom import polygon
 from odc.geo.gridspec import GridSpec
 from odc.geo.testutils import SAMPLE_WKT_WITHOUT_AUTHORITY
@@ -18,6 +18,36 @@
 # pylint: disable=comparison-with-itself,unnecessary-comprehension
 
 
+def test_gridspec_small():
+    print("Starting test for GridSpec")
+    WGS84GRID30 = GridSpec(
+        "EPSG:4326", tile_shape=(5000, 5000), resolution=0.0003, origin=XY(-180, -90)
+    )
+
+    assert WGS84GRID30.tile_shape == (5000, 5000)
+    assert WGS84GRID30.tile_size == XY(1.5, 1.5)
+
+    # Tile is at (-180 + 50*1.5) and (-90 + 50*1.5)
+    tile = (50, 50)
+    geobox = WGS84GRID30.tile_geobox(tile)
+    affine = geobox.affine
+
+    # Affine should be like this: (0.0003, 0, -105.0, 0, -0.0003, -13.5)
+    assert affine.a == 0.0003
+    assert affine.c == -105.0  # -180 + 50 * 1.5 * 0.0003
+    assert affine.f == -13.50  # -90 + 50 * 1.5 * 0.0003
+
+    # Tile is at (-180 + 200*1.5) and (-90 + 75*1.5)
+    tile = (200, 75)
+    geobox = WGS84GRID30.tile_geobox(tile)
+    affine = geobox.affine
+
+    # Affine should be like this: (0.0003, 0, 120.0, 0, -0.0003, 24.0)
+    assert affine.a == 0.0003
+    assert affine.c == 120.0  # -180 + 200 * 1.5 * 0.0003
+    assert affine.f == 24.0  # -90 + 75 * 1.5 * 0.0003
+
+
 def test_gridspec():
     gs = GridSpec(
         crs=CRS("EPSG:4326"),