remove geographic_crs and update info model

cogeotiff · Oct 8, 2024 · f0190a7 · f0190a7
1 parent bb33ffb
commit f0190a7
Show file tree

Hide file tree

Showing 8 changed files with 84 additions and 92 deletions.
diff --git a/rio_tiler/io/base.py b/rio_tiler/io/base.py
@@ -15,7 +15,7 @@
 from rasterio.rio.overview import get_maximum_overview_level
 from rasterio.warp import calculate_default_transform, transform_bounds
 
-from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
+from rio_tiler.constants import WEB_MERCATOR_TMS
 from rio_tiler.errors import (
     AssetAsBandError,
     ExpressionMixingWarning,
@@ -27,7 +27,7 @@
 from rio_tiler.models import BandStatistics, ImageData, Info, PointData
 from rio_tiler.tasks import multi_arrays, multi_points, multi_values
 from rio_tiler.types import AssetInfo, BBox, Indexes
-from rio_tiler.utils import cast_to_sequence, normalize_bounds
+from rio_tiler.utils import CRS_to_uri, cast_to_sequence, normalize_bounds
 
 
 @attr.s
@@ -48,12 +48,9 @@ class SpatialMixin:
     height: Optional[int] = attr.ib(default=None, init=False)
     width: Optional[int] = attr.ib(default=None, init=False)
 
-    geographic_crs: CRS = attr.ib(init=False, default=WGS84_CRS)
-
-    @cached_property
-    def geographic_bounds(self) -> BBox:
-        """Return dataset bounds in geographic_crs."""
-        if self.crs == self.geographic_crs:
+    def get_geographic_bounds(self, crs: CRS) -> BBox:
+        """Return Geographic Bounds for a Geographic CRS."""
+        if self.crs == crs:
             if self.bounds[1] > self.bounds[3]:
                 warnings.warn(
                     "BoundingBox of the dataset is inverted (minLat > maxLat).",
@@ -69,12 +66,7 @@ def geographic_bounds(self) -> BBox:
             return self.bounds
 
         try:
-            bounds = transform_bounds(
-                self.crs,
-                self.geographic_crs,
-                *self.bounds,
-                densify_pts=21,
-            )
+            bounds = transform_bounds(self.crs, crs, *self.bounds, densify_pts=21)
         except:  # noqa
             warnings.warn(
                 "Cannot determine bounds in geographic CRS, will default to (-180.0, -90.0, 180.0, 90.0).",
@@ -1085,9 +1077,8 @@ def _reader(band: str, **kwargs: Any) -> Info:
         bands_metadata = multi_values(bands, _reader, **kwargs)
 
         meta = {
-            "bounds": self.geographic_bounds,
-            "minzoom": self.minzoom,
-            "maxzoom": self.maxzoom,
+            "bounds": self.bounds,
+            "crs": CRS_to_uri(self.crs) or self.crs.to_wkt(),
         }
 
         # We only keep the value for the first band.

diff --git a/rio_tiler/io/rasterio.py b/rio_tiler/io/rasterio.py
@@ -35,7 +35,12 @@
 from rio_tiler.io.base import BaseReader
 from rio_tiler.models import BandStatistics, ImageData, Info, PointData
 from rio_tiler.types import BBox, Indexes, NumType, RIOResampling
-from rio_tiler.utils import _validate_shape_input, has_alpha_band, has_mask_band
+from rio_tiler.utils import (
+    CRS_to_uri,
+    _validate_shape_input,
+    has_alpha_band,
+    has_mask_band,
+)
 
 
 @attr.s
@@ -46,7 +51,6 @@ class Reader(BaseReader):
         input (str): dataset path.
         dataset (rasterio.io.DatasetReader or rasterio.io.DatasetWriter or rasterio.vrt.WarpedVRT, optional): Rasterio dataset.
         tms (morecantile.TileMatrixSet, optional): TileMatrixSet grid definition. Defaults to `WebMercatorQuad`.
-        geographic_crs (rasterio.crs.CRS, optional): CRS to use as geographic coordinate system. Defaults to WGS84.
         colormap (dict, optional): Overwrite internal colormap.
         options (dict, optional): Options to forward to low-level reader methods.
 
@@ -75,7 +79,6 @@ class Reader(BaseReader):
     )
 
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
-    geographic_crs: CRS = attr.ib(default=WGS84_CRS)
 
     colormap: Dict = attr.ib(default=None)
 
@@ -177,9 +180,8 @@ def _get_descr(ix):
             nodata_type = "None"
 
         meta = {
-            "bounds": self.geographic_bounds,
-            "minzoom": self.minzoom,
-            "maxzoom": self.maxzoom,
+            "bounds": self.bounds,
+            "crs": CRS_to_uri(self.crs) or self.crs.to_wkt(),
             "band_metadata": [
                 (f"b{ix}", self.dataset.tags(ix)) for ix in self.dataset.indexes
             ],
@@ -596,8 +598,6 @@ class ImageReader(Reader):
     tms: TileMatrixSet = attr.ib(init=False)
 
     crs: CRS = attr.ib(init=False, default=None)
-    geographic_crs: CRS = attr.ib(init=False, default=None)
-
     transform: Affine = attr.ib(init=False)
 
     def __attrs_post_init__(self):

diff --git a/rio_tiler/io/stac.py b/rio_tiler/io/stac.py
@@ -13,7 +13,6 @@
 from cachetools import LRUCache, cached
 from cachetools.keys import hashkey
 from morecantile import TileMatrixSet
-from rasterio.crs import CRS
 from rasterio.transform import array_bounds
 
 from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
@@ -198,7 +197,6 @@ class STACReader(MultiBaseReader):
         tms (morecantile.TileMatrixSet, optional): TileMatrixSet grid definition. Defaults to `WebMercatorQuad`.
         minzoom (int, optional): Set minzoom for the tiles.
         maxzoom (int, optional): Set maxzoom for the tiles.
-        geographic_crs (rasterio.crs.CRS, optional): CRS to use as geographic coordinate system. Defaults to WGS84.
         include_assets (set of string, optional): Only Include specific assets.
         exclude_assets (set of string, optional): Exclude specific assets.
         include_asset_types (set of string, optional): Only include some assets base on their type.
@@ -234,8 +232,6 @@ class STACReader(MultiBaseReader):
     minzoom: int = attr.ib(default=None)
     maxzoom: int = attr.ib(default=None)
 
-    geographic_crs: CRS = attr.ib(default=WGS84_CRS)
-
     include_assets: Optional[Set[str]] = attr.ib(default=None)
     exclude_assets: Optional[Set[str]] = attr.ib(default=None)
 

diff --git a/rio_tiler/io/xarray.py b/rio_tiler/io/xarray.py
@@ -23,7 +23,7 @@
 from rio_tiler.io.base import BaseReader
 from rio_tiler.models import BandStatistics, ImageData, Info, PointData
 from rio_tiler.types import BBox, NoData, WarpResampling
-from rio_tiler.utils import _validate_shape_input
+from rio_tiler.utils import CRS_to_uri, _validate_shape_input
 
 try:
     import xarray
@@ -43,7 +43,6 @@ class XarrayReader(BaseReader):
     Attributes:
         dataset (xarray.DataArray): Xarray DataArray dataset.
         tms (morecantile.TileMatrixSet, optional): TileMatrixSet grid definition. Defaults to `WebMercatorQuad`.
-        geographic_crs (rasterio.crs.CRS, optional): CRS to use as geographic coordinate system. Defaults to WGS84.
 
     Examples:
         >>> ds = xarray.open_dataset(
@@ -62,7 +61,6 @@ class XarrayReader(BaseReader):
     input: xarray.DataArray = attr.ib()
 
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
-    geographic_crs: CRS = attr.ib(default=WGS84_CRS)
 
     _dims: List = attr.ib(init=False, factory=list)
 
@@ -120,9 +118,8 @@ def info(self) -> Info:
         metadata = [band.attrs for d in self._dims for band in self.input[d]]
 
         meta = {
-            "bounds": self.geographic_bounds,
-            "minzoom": self.minzoom,
-            "maxzoom": self.maxzoom,
+            "bounds": self.bounds,
+            "crs": CRS_to_uri(self.crs),
             "band_metadata": [(f"b{ix}", v) for ix, v in enumerate(metadata, 1)],
             "band_descriptions": [(f"b{ix}", v) for ix, v in enumerate(bands, 1)],
             "dtype": str(self.input.dtype),

diff --git a/rio_tiler/models.py b/rio_tiler/models.py
@@ -61,16 +61,10 @@ class Bounds(RioTilerBaseModel):
     """Dataset Bounding box"""
 
     bounds: BoundingBox
+    crs: str
 
 
-class SpatialInfo(Bounds):
-    """Dataset SpatialInfo"""
-
-    minzoom: int
-    maxzoom: int
-
-
-class Info(SpatialInfo):
+class Info(Bounds):
     """Dataset Info."""
 
     band_metadata: List[Tuple[str, Dict]]

diff --git a/rio_tiler/utils.py b/rio_tiler/utils.py
@@ -789,3 +789,21 @@ def cast_to_sequence(val: Optional[Any] = None) -> Sequence:
         val = (val,)
 
     return val
+
+
+def CRS_to_uri(crs: CRS) -> Optional[str]:
+    """Convert CRS to URI.
+
+    Code adapted from https://github.com/developmentseed/morecantile/blob/1829fe12408e4a1feee7493308f3f02257ef4caf/morecantile/models.py#L148-L161
+    """
+    # attempt to grab the authority, version, and code from the CRS
+    if authority_code := crs.to_authority(confidence_threshold=70):
+        version = "0"
+        authority, code = authority_code
+        # if we have a version number in the authority, split it out
+        if "_" in authority:
+            authority, version = authority.split("_")
+
+        return f"http://www.opengis.net/def/crs/{authority}/{version}/{code}"
+
+    return None
diff --git a/tests/test_io_rasterio.py b/tests/test_io_rasterio.py
@@ -1,18 +1,17 @@
 """tests rio_tiler.io.rasterio.Reader"""
 
 import os
-import warnings
 from io import BytesIO
 from typing import Any, Dict
 
-import attr
 import morecantile
 import numpy
 import pytest
 import rasterio
 from morecantile import TileMatrixSet
 from pyproj import CRS
 from rasterio import transform
+from rasterio.crs import CRS as rioCRS
 from rasterio.features import bounds as featureBounds
 from rasterio.io import MemoryFile
 from rasterio.vrt import WarpedVRT
@@ -82,6 +81,10 @@ def test_info_valid():
     """Should work as expected (get file info)"""
     with Reader(COG_SCALE) as src:
         meta = src.info()
+        assert meta.bounds == src.bounds == src.dataset.bounds
+        crs = meta.crs
+        assert rioCRS.from_user_input(crs) == src.crs
+
         assert meta.scales
         assert meta.offsets
         assert not meta.colormap
@@ -92,6 +95,11 @@ def test_info_valid():
         assert meta.driver
 
     with Reader(COG_CMAP) as src:
+        meta = src.info()
+        assert meta.bounds == src.bounds == src.dataset.bounds
+        crs = meta.crs
+        assert rioCRS.from_user_input(crs) == src.crs
+
         assert src.colormap
         meta = src.info()
         assert meta["colormap"]
@@ -106,8 +114,7 @@ def test_info_valid():
     with Reader(COG_TAGS) as src:
         meta = src.info()
         assert meta.bounds
-        assert meta.minzoom
-        assert meta.maxzoom
+        assert meta.crs
         assert meta.band_descriptions
         assert meta.dtype == "int16"
         assert meta.colorinterp == ["gray"]
@@ -894,26 +901,30 @@ def test_nonearthbody():
             assert src.maxzoom == 24
 
     # Warns because of zoom level in WebMercator can't be defined
-    with pytest.warns(UserWarning) as w:
-        with Reader(COG_EUROPA, geographic_crs=EUROPA_SPHERE) as src:
-            assert src.info()
-            assert len(w) == 2
+    with Reader(COG_EUROPA) as src:
+        assert src.info()
+        meta = src.info()
+        assert meta.bounds == src.bounds == src.dataset.bounds
+        crs = meta.crs
+        assert rioCRS.from_user_input(crs) == src.crs
 
-            img = src.read()
-            assert numpy.array_equal(img.data, src.dataset.read(indexes=(1,)))
-            assert img.width == src.dataset.width
-            assert img.height == src.dataset.height
-            assert img.count == src.dataset.count
+        assert src.get_geographic_bounds(EUROPA_SPHERE)
+
+        img = src.read()
+        assert numpy.array_equal(img.data, src.dataset.read(indexes=(1,)))
+        assert img.width == src.dataset.width
+        assert img.height == src.dataset.height
+        assert img.count == src.dataset.count
 
-            img = src.preview()
-            assert img.bounds == src.bounds
+        img = src.preview()
+        assert img.bounds == src.bounds
 
-            part = src.part(src.bounds, bounds_crs=src.crs)
-            assert part.bounds == src.bounds
+        part = src.part(src.bounds, bounds_crs=src.crs)
+        assert part.bounds == src.bounds
 
-            lon = (src.bounds[0] + src.bounds[2]) / 2
-            lat = (src.bounds[1] + src.bounds[3]) / 2
-            assert src.point(lon, lat, coord_crs=src.crs).data[0] is not None
+        lon = (src.bounds[0] + src.bounds[2]) / 2
+        lat = (src.bounds[1] + src.bounds[3]) / 2
+        assert src.point(lon, lat, coord_crs=src.crs).data[0] is not None
 
     with pytest.warns(UserWarning):
         europa_crs = CRS.from_authority("ESRI", 104915)
@@ -923,7 +934,7 @@ def test_nonearthbody():
             matrix_scale=[2, 1],
         )
 
-    with Reader(COG_EUROPA, tms=tms, geographic_crs=EUROPA_SPHERE) as src:
+    with Reader(COG_EUROPA, tms=tms) as src:
         assert src.info()
         assert src.minzoom == 4
         assert src.maxzoom == 6
@@ -958,28 +969,8 @@ def test_nonearth_custom():
         geographic_crs=MARS2000_SPHERE,
     )
 
-    @attr.s
-    class MarsReader(Reader):
-        """Use custom geographic CRS."""
-
-        geographic_crs: rasterio.crs.CRS = attr.ib(
-            init=False,
-            default=rasterio.crs.CRS.from_proj4("+proj=longlat +R=3396190 +no_defs"),
-        )
-
-    with warnings.catch_warnings():
-        with MarsReader(COG_MARS, tms=mars_tms) as src:
-            assert src.geographic_bounds[0] > -180
-
-    with warnings.catch_warnings():
-        with Reader(
-            COG_MARS,
-            tms=mars_tms,
-            geographic_crs=rasterio.crs.CRS.from_proj4(
-                "+proj=longlat +R=3396190 +no_defs"
-            ),
-        ) as src:
-            assert src.geographic_bounds[0] > -180
+    with Reader(COG_MARS, tms=mars_tms) as src:
+        assert src.get_geographic_bounds(MARS2000_SPHERE)[0] > -180
 
 
 def test_tms_tilesize_and_zoom():
@@ -1121,7 +1112,10 @@ def test_inverted_latitude():
     """Test working with inverted Latitude."""
     with pytest.warns(UserWarning):
         with Reader(COG_INVERTED) as src:
-            assert src.geographic_bounds[1] < src.geographic_bounds[3]
+            assert (
+                src.get_geographic_bounds(WGS84_CRS)[1]
+                < src.get_geographic_bounds(WGS84_CRS)[3]
+            )
 
     with pytest.warns(UserWarning):
         with Reader(COG_INVERTED) as src:

diff --git a/tests/test_io_xarray.py b/tests/test_io_xarray.py
@@ -7,6 +7,7 @@
 import pytest
 import rioxarray
 import xarray
+from rasterio.crs import CRS as rioCRS
 
 from rio_tiler.errors import InvalidGeographicBounds, MissingCRS
 from rio_tiler.io import XarrayReader
@@ -28,9 +29,11 @@ def test_xarray_reader():
 
     data.rio.write_crs("epsg:4326", inplace=True)
     with XarrayReader(data) as dst:
+        assert dst.minzoom == dst.maxzoom == 0
         info = dst.info()
-        assert info.minzoom == 0
-        assert info.maxzoom == 0
+        assert info.bounds == dst.bounds
+        crs = info.crs
+        assert rioCRS.from_user_input(crs) == dst.crs
         assert info.band_metadata == [("b1", {})]
         assert info.band_descriptions == [("b1", "2022-01-01T00:00:00.000000000")]
         assert info.height == 33
@@ -123,9 +126,8 @@ def test_xarray_reader():
 
     data.rio.write_crs("epsg:4326", inplace=True)
     with XarrayReader(data) as dst:
-        info = dst.info()
-        assert info.minzoom == 5
-        assert info.maxzoom == 7
+        assert dst.minzoom == 5
+        assert dst.maxzoom == 7
 
 
 def test_xarray_reader_external_nodata():