refactor SpatialMixin.tile_exists to avoid coordinates overflow

CHANGES.md

@@ -1,4 +1,8 @@
 
+# 3.0.0a4 (2021-11-10)
+
+* refactor `SpatialMixin.tile_exists` to compare the bounds in the dataset's coordinate system to avoid coordinates overflow (a TMS CRS bounds can be smaller than the dataset CRS bounds) (https://github.com/cogeotiff/rio-tiler/pull/455)
+
 # 3.0.0a3 (2021-11-03)
 
 * Reader's `info` and `statistics` methods to default to available `bands` or `assets` if not provided (https://github.com/cogeotiff/rio-tiler/pull/451)

rio_tiler/io/base.py

@@ -6,6 +6,7 @@
 from typing import Any, Coroutine, Dict, List, Optional, Sequence, Tuple, Type, Union
 
 import attr
+import numpy
 from morecantile import Tile, TileMatrixSet
 from rasterio.crs import CRS
 from rasterio.warp import transform_bounds
@@ -74,31 +75,37 @@ def tile_exists(self, tile_x: int, tile_y: int, tile_z: int) -> bool:
             bool: True if the tile intersects the dataset bounds.
 
         """
+        # bounds in TileMatrixSet's CRS
         tile_bounds = self.tms.xy_bounds(Tile(x=tile_x, y=tile_y, z=tile_z))
 
+        # Transform the bounds to the dataset's CRS
         try:
-            dataset_bounds = transform_bounds(
-                self.crs,
+            tile_bounds = transform_bounds(
                 self.tms.rasterio_crs,
-                *self.bounds,
+                self.crs,
+                *tile_bounds,
                 densify_pts=21,
             )
         except:  # noqa
             # HACK: gdal will first throw an error for invalid transformation
             # but if retried it will then pass.
             # Note: It might return `+/-inf` values
-            dataset_bounds = transform_bounds(
-                self.crs,
+            tile_bounds = transform_bounds(
                 self.tms.rasterio_crs,
-                *self.bounds,
+                self.crs,
+                *tile_bounds,
                 densify_pts=21,
             )
 
+        # If tile_bounds has non-finite value in the dataset CRS we return True
+        if not all(numpy.isfinite(tile_bounds)):
+            return True
+
         return (
-            (tile_bounds[0] < dataset_bounds[2])
-            and (tile_bounds[2] > dataset_bounds[0])
-            and (tile_bounds[3] > dataset_bounds[1])
-            and (tile_bounds[1] < dataset_bounds[3])
+            (tile_bounds[0] < self.bounds[2])
+            and (tile_bounds[2] > self.bounds[0])
+            and (tile_bounds[3] > self.bounds[1])
+            and (tile_bounds[1] < self.bounds[3])
         )
 
 

tests/test_io_cogeo.py

@@ -4,6 +4,7 @@
 from io import BytesIO
 from typing import Any, Dict
 
+import morecantile
 import numpy
 import pytest
 import rasterio
@@ -760,6 +761,12 @@ def test_fullEarth():
         img = cog.tile(127, 42, 7, tilesize=64)
         assert img.data.shape == (1, 64, 64)
 
+    with COGReader(
+        COG_EARTH, tms=morecantile.tms.get("EuropeanETRS89_LAEAQuad")
+    ) as cog:
+        img = cog.tile(0, 0, 1, tilesize=64)
+        assert img.data.shape == (1, 64, 64)
+
 
 def test_read():
     """Should read the entire dataset."""