allow wgs84 geographic crs (#458)

* allow wgs84 geographic crs * update changelog * surface geographic_crs in Readers * Update docs/v3_migration.md
cogeotiff · Nov 16, 2021 · fdb6836 · fdb6836
1 parent 6470b54
commit fdb6836
Show file tree

Hide file tree

Showing 7 changed files with 100 additions and 26 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -1,4 +1,8 @@
 
+# 3.0.0a5 (TBD)
+
+* allow the definition of `geographic_crs` used in the `geographic_bounds` property (https://github.com/cogeotiff/rio-tiler/pull/458)
+
 # 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)

diff --git a/docs/v3_migration.md b/docs/v3_migration.md
@@ -29,6 +29,34 @@ with COGReader("my_tif.tif") as cog:
     print(cog.geographic_bounds)  # bounds in WGS84 projection
 ```
 
+Note a `geographic_crs` attribute is available in `COGReader` and `STACReader` to control which CRS to use for the transformation from the dataset's CRS.  `geographic_crs` is outside the `__init__` method for Abstract Base Classes (e.g `BaseReader`)
+
+```python
+MARS2000_SPHERE = CRS.from_proj4("+proj=longlat +R=3396190 +no_defs")
+MARS_MERCATOR = CRS.from_proj4(
+    "+proj=merc +R=3396190 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs"
+)
+MARS_TMS = TileMatrixSet.custom(
+    [
+        -179.9999999999996,
+        -85.05112877980656,
+        179.9999999999996,
+        85.05112877980656,
+    ],
+    MARS_MERCATOR,
+    extent_crs=MARS2000_SPHERE,
+    title="Web Mercator Mars",
+    geographic_crs=MARS2000_SPHERE,
+)
+with COGReader(
+    "martian_dataset.tif",
+    tms=MARS_TMS,
+    geographic_crs=rasterio.crs.CRS.from_proj4("+proj=longlat +R=3396190 +no_defs"),
+) as cog:
+    ...
+```
+
+
 ## No more `max_size`
 
 We've removed the default option `max_size=1024` in `BaseReader.part` and `BaseReader.feature` to return the highest resolution dataset by default.

diff --git a/rio_tiler/io/base.py b/rio_tiler/io/base.py
@@ -38,19 +38,22 @@ class SpatialMixin:
     """
 
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
+
     minzoom: int = attr.ib(init=False)
     maxzoom: int = attr.ib(init=False)
 
     bounds: BBox = attr.ib(init=False)
     crs: CRS = attr.ib(init=False)
 
+    geographic_crs: CRS = attr.ib(init=False, default=WGS84_CRS)
+
     @property
     def geographic_bounds(self) -> BBox:
         """return bounds in WGS84."""
         try:
             bounds = transform_bounds(
                 self.crs,
-                WGS84_CRS,
+                self.geographic_crs,
                 *self.bounds,
                 densify_pts=21,
             )
@@ -116,12 +119,6 @@ class BaseReader(SpatialMixin, metaclass=abc.ABCMeta):
     input: Any = attr.ib()
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
 
-    minzoom: int = attr.ib(init=False)
-    maxzoom: int = attr.ib(init=False)
-
-    bounds: BBox = attr.ib(init=False)
-    crs: CRS = attr.ib(init=False)
-
     def __enter__(self):
         """Support using with Context Managers."""
         return self
@@ -223,12 +220,6 @@ class AsyncBaseReader(SpatialMixin, metaclass=abc.ABCMeta):
     input: Any = attr.ib()
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
 
-    minzoom: int = attr.ib(init=False)
-    maxzoom: int = attr.ib(init=False)
-
-    bounds: BBox = attr.ib(init=False)
-    crs: CRS = attr.ib(init=False)
-
     async def __aenter__(self):
         """Support using with Context Managers."""
         return self
@@ -339,6 +330,7 @@ class MultiBaseReader(BaseReader, metaclass=abc.ABCMeta):
     This Reader is suited for dataset that are composed of multiple assets (e.g. STAC).
 
     Attributes:
+        input (any): input data.
         reader (rio_tiler.io.BaseReader): reader.
         reader_options (dict, option): options to forward to the reader. Defaults to `{}`.
         tms (morecantile.TileMatrixSet, optional): TileMatrixSet grid definition. Defaults to `WebMercatorQuad`.
@@ -352,12 +344,6 @@ class MultiBaseReader(BaseReader, metaclass=abc.ABCMeta):
     reader_options: Dict = attr.ib(factory=dict)
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
 
-    minzoom: int = attr.ib(init=False)
-    maxzoom: int = attr.ib(init=False)
-
-    bounds: BBox = attr.ib(init=False)
-    crs: CRS = attr.ib(init=False)
-
     assets: Sequence[str] = attr.ib(init=False)
 
     @abc.abstractmethod
@@ -787,6 +773,7 @@ class MultiBandReader(BaseReader, metaclass=abc.ABCMeta):
     This Reader is suited for dataset that stores spectral bands as separate files  (e.g. Sentinel 2).
 
     Attributes:
+        input (any): input data.
         reader (rio_tiler.io.BaseReader): reader.
         reader_options (dict, option): options to forward to the reader. Defaults to `{}`.
         tms (morecantile.TileMatrixSet, optional): TileMatrixSet grid definition. Defaults to `WebMercatorQuad`.
@@ -800,12 +787,6 @@ class MultiBandReader(BaseReader, metaclass=abc.ABCMeta):
     reader_options: Dict = attr.ib(factory=dict)
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
 
-    minzoom: int = attr.ib(init=False)
-    maxzoom: int = attr.ib(init=False)
-
-    bounds: BBox = attr.ib(init=False)
-    crs: CRS = attr.ib(init=False)
-
     bands: Sequence[str] = attr.ib(init=False)
 
     @abc.abstractmethod

diff --git a/rio_tiler/io/cogeo.py b/rio_tiler/io/cogeo.py
@@ -44,10 +44,11 @@ class COGReader(BaseReader):
         input (str): Cloud Optimized GeoTIFF path.
         dataset (rasterio.io.DatasetReader or rasterio.io.DatasetWriter or rasterio.vrt.WarpedVRT, optional): Rasterio dataset.
         bounds (tuple): Dataset bounds (left, bottom, right, top).
-        crs (rasterio.crs.CRS): Dataset crs.
+        crs (rasterio.crs.CRS): Dataset CRS.
         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.
         colormap (dict, optional): Overwrite internal colormap.
         nodata (int or float or str, optional): Global options, overwrite internal nodata value.
         unscale (bool, optional): Global options, apply internal scale and offset on all read operations.
@@ -83,6 +84,8 @@ class COGReader(BaseReader):
     minzoom: int = attr.ib(default=None)
     maxzoom: int = attr.ib(default=None)
 
+    geographic_crs: CRS = attr.ib(default=WGS84_CRS)
+
     colormap: Dict = attr.ib(default=None)
 
     # Define global options to be forwarded to functions reading the data (e.g `rio_tiler.reader.read`)
@@ -672,9 +675,13 @@ class GCPCOGReader(COGReader):
     src_dataset: Union[DatasetReader, DatasetWriter, MemoryFile, WarpedVRT] = attr.ib(
         default=None
     )
+
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
     minzoom: int = attr.ib(default=None)
     maxzoom: int = attr.ib(default=None)
+
+    geographic_crs: CRS = attr.ib(default=WGS84_CRS)
+
     colormap: Dict = attr.ib(default=None)
 
     # Define global options to be forwarded to functions reading the data (e.g `rio_tiler.reader.read`)

diff --git a/rio_tiler/io/stac.py b/rio_tiler/io/stac.py
@@ -9,6 +9,7 @@
 import httpx
 import pystac
 from morecantile import TileMatrixSet
+from rasterio.crs import CRS
 
 from ..constants import WEB_MERCATOR_TMS, WGS84_CRS
 from ..errors import InvalidAssetName, MissingAssets
@@ -126,6 +127,7 @@ class STACReader(MultiBaseReader):
         item (dict or pystac.Item, STAC): Stac Item.
         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 (set of string, optional): Only Include specific assets.
         exclude (set of string, optional): Exclude specific assets.
         include_asset_types (set of string, optional): Only include some assets base on their type.
@@ -160,6 +162,8 @@ 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/tests/fixtures/cog_hirise_mars.tif b/tests/fixtures/cog_hirise_mars.tif
diff --git a/tests/test_io_cogeo.py b/tests/test_io_cogeo.py
@@ -4,6 +4,7 @@
 from io import BytesIO
 from typing import Any, Dict
 
+import attr
 import morecantile
 import numpy
 import pytest
@@ -38,6 +39,7 @@
 COG_EARTH = os.path.join(PREFIX, "cog_fullearth.tif")
 GEOTIFF = os.path.join(PREFIX, "nocog.tif")
 COG_EUROPA = os.path.join(PREFIX, "cog_nonearth.tif")
+COG_MARS = os.path.join(PREFIX, "cog_hirise_mars.tif")
 
 KEY_ALPHA = "hro_sources/colorado/201404_13SED190110_201404_0x1500m_CL_1_alpha.tif"
 COG_ALPHA = os.path.join(PREFIX, "my-bucket", KEY_ALPHA)
@@ -849,3 +851,51 @@ def test_nonearthbody():
             assert img.crs == tms.rasterio_crs
 
             assert len(warnings) == 0
+
+
+def test_nonearth_custom():
+    """Test Custom geographic_crs."""
+    MARS2000_SPHERE = CRS.from_proj4("+proj=longlat +R=3396190 +no_defs")
+    MARS_MERCATOR = CRS.from_proj4(
+        "+proj=merc +R=3396190 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs"
+    )
+
+    mars_tms = TileMatrixSet.custom(
+        [
+            -179.9999999999996,
+            -85.05112877980656,
+            179.9999999999996,
+            85.05112877980656,
+        ],
+        MARS_MERCATOR,
+        extent_crs=MARS2000_SPHERE,
+        title="Web Mercator Mars",
+        geographic_crs=MARS2000_SPHERE,
+    )
+
+    @attr.s
+    class MarsReader(COGReader):
+        """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 pytest.warns(None) as warnings:
+        with MarsReader(COG_MARS, tms=mars_tms) as cog:
+            assert cog.geographic_bounds[0] > -180
+
+    assert len(warnings) == 0
+
+    with pytest.warns(None) as warnings:
+        with COGReader(
+            COG_MARS,
+            tms=mars_tms,
+            geographic_crs=rasterio.crs.CRS.from_proj4(
+                "+proj=longlat +R=3396190 +no_defs"
+            ),
+        ) as cog:
+            assert cog.geographic_bounds[0] > -180
+
+    assert len(warnings) == 0