add warning when non wgs84 geographic crs and optional wgs84 default

developmentseed · Nov 8, 2024 · 31994e7 · 31994e7
1 parent 0ab27a7
commit 31994e7
Show file tree

Hide file tree

Showing 6 changed files with 92 additions and 47 deletions.
diff --git a/.gitignore b/.gitignore
@@ -134,3 +134,5 @@ dmypy.json
 # VSCode
 .vscode
 .vscode/
+
+.notebooks/
diff --git a/morecantile/errors.py b/morecantile/errors.py
@@ -35,3 +35,7 @@ class InvalidZoomError(MorecantileError):
 
 class DeprecationError(MorecantileError):
     """Raised when TMS version is not 2.0"""
+
+
+class NonWGS84GeographicCRS(UserWarning):
+    """Geographic CRS is not EPSG:4326."""
diff --git a/morecantile/models.py b/morecantile/models.py
@@ -1,6 +1,7 @@
 """Pydantic modules for OGC TileMatrixSets (https://www.ogc.org/standards/tms)"""
 
 import math
+import os
 import warnings
 from functools import cached_property
 from typing import Any, Dict, Iterator, List, Literal, Optional, Sequence, Tuple, Union
@@ -23,6 +24,7 @@
 from morecantile.errors import (
     DeprecationError,
     InvalidZoomError,
+    NonWGS84GeographicCRS,
     NoQuadkeySupport,
     PointOutsideTMSBounds,
     QuadKeyError,
@@ -41,6 +43,14 @@
 BoundsType = Tuple[NumType, NumType]
 LL_EPSILON = 1e-11
 axesInfo = Annotated[List[str], Field(min_length=2, max_length=2)]
+WGS84_CRS = pyproj.CRS.from_epsg(4326)
+
+IGNORE_NON_WGS84_GEOGRAPHIC = os.environ.get(
+    "MORECANTILE_IGNORE_NON_WGS84_GEOGRAPHIC", "false"
+).lower() in ["true", "on", "1", "yes"]
+DEFAULT_WGS84_GEOGRAPHIC = os.environ.get(
+    "MORECANTILE_WGS84_GEOGRAPHIC", "false"
+).lower() in ["true", "on", "1", "yes"]
 
 
 class CRSUri(BaseModel):
@@ -485,6 +495,7 @@ class TileMatrixSet(BaseModel, arbitrary_types_allowed=True):
     ]
 
     # Private attributes
+    _geographic_crs: pyproj.CRS = PrivateAttr()
     _to_geographic: pyproj.Transformer = PrivateAttr()
     _from_geographic: pyproj.Transformer = PrivateAttr()
 
@@ -499,11 +510,22 @@ def __init__(self, **data):
         }
 
         try:
+            self._geographic_crs = (
+                WGS84_CRS
+                if DEFAULT_WGS84_GEOGRAPHIC
+                else self.crs._pyproj_crs.geodetic_crs
+            )
+            if not IGNORE_NON_WGS84_GEOGRAPHIC and self._geographic_crs != WGS84_CRS:
+                warnings.warn(
+                    f"`{self.id}` TMS's CRS doesn't use EPSG:4326 as geographic CRS but {self._geographic_crs}",
+                    NonWGS84GeographicCRS,
+                )
+
             self._to_geographic = pyproj.Transformer.from_crs(
-                self.crs._pyproj_crs, self.crs._pyproj_crs.geodetic_crs, always_xy=True
+                self.crs._pyproj_crs, self._geographic_crs, always_xy=True
             )
             self._from_geographic = pyproj.Transformer.from_crs(
-                self.crs._pyproj_crs.geodetic_crs, self.crs._pyproj_crs, always_xy=True
+                self._geographic_crs, self.crs._pyproj_crs, always_xy=True
             )
         except ProjError:
             warnings.warn(
@@ -555,7 +577,7 @@ def __repr__(self):
     @cached_property
     def geographic_crs(self) -> pyproj.CRS:
         """Return the TMS's geographic CRS."""
-        return self.crs._pyproj_crs.geodetic_crs
+        return self._geographic_crs
 
     @cached_property
     def rasterio_crs(self):
@@ -565,7 +587,7 @@ def rasterio_crs(self):
     @cached_property
     def rasterio_geographic_crs(self):
         """Return the geographic CRS as a rasterio CRS."""
-        return to_rasterio_crs(self.crs._pyproj_crs.geodetic_crs)
+        return to_rasterio_crs(self._geographic_crs)
 
     @property
     def minzoom(self) -> int:

diff --git a/pyproject.toml b/pyproject.toml
@@ -117,6 +117,12 @@ ignore = [
 [tool.pytest.ini_options]
 filterwarnings = [
     "ignore:You will likely lose important projection*:UserWarning",
+    "ignore:`CanadianNAD83_LCC`*:morecantile.errors.NonWGS84GeographicCRS",
+    "ignore:`EuropeanETRS89_LAEAQuad`*:morecantile.errors.NonWGS84GeographicCRS",
+    "ignore:`WorldCRS84Quad`*:morecantile.errors.NonWGS84GeographicCRS",
+    "ignore:`NZTM2000Quad`*:morecantile.errors.NonWGS84GeographicCRS",
+    "ignore:`LINZAntarticaMapTilegrid`*:morecantile.errors.NonWGS84GeographicCRS",
+    "ignore::morecantile.errors.PointOutsideTMSBounds",
 ]
 
 [tool.bumpversion]

diff --git a/tests/test_models.py b/tests/test_models.py
@@ -12,7 +12,7 @@
 
 import morecantile
 from morecantile.commons import Tile
-from morecantile.errors import InvalidIdentifier
+from morecantile.errors import InvalidIdentifier, NonWGS84GeographicCRS
 from morecantile.models import CRS, CRSWKT, CRSUri, TileMatrix, TileMatrixSet
 
 data_dir = os.path.join(os.path.dirname(__file__), "../morecantile/data")
@@ -214,11 +214,12 @@ def test_custom_tms_decimation():
     extent = (238170, 4334121, 377264, 4473215)
     left, bottom, right, top = extent
     for decimation_base in [2, 3, 4, 5]:
-        custom_tms = TileMatrixSet.custom(
-            extent,
-            pyproj.CRS.from_epsg(6342),
-            decimation_base=decimation_base,
-        )
+        with pytest.warns(NonWGS84GeographicCRS):
+            custom_tms = TileMatrixSet.custom(
+                extent,
+                pyproj.CRS.from_epsg(6342),
+                decimation_base=decimation_base,
+            )
 
         if decimation_base == 2:
             assert custom_tms.is_quadtree
@@ -311,7 +312,8 @@ def test_schema():
         "+proj=stere +lat_0=90 +lon_0=0 +k=2 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs"
     )
     extent = [-13584760.000, -13585240.000, 13585240.000, 13584760.000]
-    tms = morecantile.TileMatrixSet.custom(extent, crs, id="MarsNPolek2MOLA5k")
+    with pytest.warns(NonWGS84GeographicCRS):
+        tms = morecantile.TileMatrixSet.custom(extent, crs, id="MarsNPolek2MOLA5k")
     assert tms.model_json_schema()
     assert tms.model_dump(exclude_none=True)
     json_doc = json.loads(tms.model_dump_json(exclude_none=True))
@@ -336,18 +338,19 @@ def test_mars_tms():
         "+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"
     )
 
-    # same boundaries as Earth mercator
-    mars_tms = TileMatrixSet.custom(
-        [
-            -179.9999999999996,
-            -85.05112877980656,
-            179.9999999999996,
-            85.05112877980656,
-        ],
-        MARS_MERCATOR,
-        extent_crs=MARS2000_SPHERE,
-        title="Web Mercator Mars",
-    )
+    with pytest.warns(NonWGS84GeographicCRS):
+        # same boundaries as Earth mercator
+        mars_tms = TileMatrixSet.custom(
+            [
+                -179.9999999999996,
+                -85.05112877980656,
+                179.9999999999996,
+                85.05112877980656,
+            ],
+            MARS_MERCATOR,
+            extent_crs=MARS2000_SPHERE,
+            title="Web Mercator Mars",
+        )
     assert mars_tms.geographic_crs == MARS2000_SPHERE
 
     pos = (35, 40, 3)
@@ -374,12 +377,15 @@ def test_mars_local_tms():
     SYRTIS_TM = pyproj.CRS.from_proj4(
         "+proj=tmerc +lat_0=17 +lon_0=76.5 +k=0.9996 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs"
     )
-    # 100km grid centered on 17N, 76.5E
-    syrtis_tms = TileMatrixSet.custom(
-        [-5e5, -5e5, 5e5, 5e5],
-        SYRTIS_TM,
-        title="Web Mercator Mars",
-    )
+
+    with pytest.warns(NonWGS84GeographicCRS):
+
+        # 100km grid centered on 17N, 76.5E
+        syrtis_tms = TileMatrixSet.custom(
+            [-5e5, -5e5, 5e5, 5e5],
+            SYRTIS_TM,
+            title="Web Mercator Mars",
+        )
     assert SYRTIS_TM == syrtis_tms.crs._pyproj_crs
     assert syrtis_tms.geographic_crs
     assert syrtis_tms.model_dump(mode="json")
@@ -398,12 +404,13 @@ def test_mars_local_tms():
 def test_mars_tms_construction():
     mars_sphere_crs = pyproj.CRS.from_user_input("IAU_2015:49900")
     extent = [-180.0, -90.0, 180.0, 90.0]
-    mars_tms = morecantile.TileMatrixSet.custom(
-        extent,
-        crs=mars_sphere_crs,
-        id="MarsGeographicCRS",
-        matrix_scale=[2, 1],
-    )
+    with pytest.warns(NonWGS84GeographicCRS):
+        mars_tms = morecantile.TileMatrixSet.custom(
+            extent,
+            crs=mars_sphere_crs,
+            id="MarsGeographicCRS",
+            matrix_scale=[2, 1],
+        )
     assert "4326" not in mars_tms.geographic_crs.to_wkt()
     assert "4326" not in mars_tms.rasterio_geographic_crs.to_wkt()
     assert mars_tms.xy_bbox.left == pytest.approx(-180.0)
@@ -421,11 +428,12 @@ def test_mars_web_mercator_long_lat():
         10669445.554195119,
         10669445.554195119,
     ]
-    mars_tms_wm = morecantile.TileMatrixSet.custom(
-        extent_wm,
-        crs=crs_mars_web_mercator,
-        id="MarsWebMercator",
-    )
+    with pytest.warns(NonWGS84GeographicCRS):
+        mars_tms_wm = morecantile.TileMatrixSet.custom(
+            extent_wm,
+            crs=crs_mars_web_mercator,
+            id="MarsWebMercator",
+        )
     assert "4326" not in mars_tms_wm.geographic_crs.to_wkt()
     assert "4326" not in mars_tms_wm.rasterio_geographic_crs.to_wkt()
     assert mars_tms_wm.bbox.left == pytest.approx(-180.0)
@@ -439,12 +447,13 @@ def test_mars_web_mercator_long_lat():
         85.05112877980656,
     ]
     mars_sphere_crs = pyproj.CRS.from_user_input("IAU_2015:49900")
-    mars_tms_wm_geog_ext = morecantile.TileMatrixSet.custom(
-        extent_wm_geog,
-        extent_crs=mars_sphere_crs,
-        crs=crs_mars_web_mercator,
-        id="MarsWebMercator",
-    )
+    with pytest.warns(NonWGS84GeographicCRS):
+        mars_tms_wm_geog_ext = morecantile.TileMatrixSet.custom(
+            extent_wm_geog,
+            extent_crs=mars_sphere_crs,
+            crs=crs_mars_web_mercator,
+            id="MarsWebMercator",
+        )
     assert mars_tms_wm_geog_ext.bbox.left == pytest.approx(-180.0)
     assert mars_tms_wm_geog_ext.bbox.bottom == pytest.approx(-85.0511287)
     assert mars_tms_wm_geog_ext.bbox.right == pytest.approx(180.0)

diff --git a/tests/test_morecantile.py b/tests/test_morecantile.py
@@ -11,6 +11,7 @@
 from morecantile.errors import (
     InvalidIdentifier,
     InvalidZoomError,
+    NonWGS84GeographicCRS,
     PointOutsideTMSBounds,
 )
 from morecantile.utils import is_power_of_two, meters_per_unit
@@ -439,7 +440,8 @@ def test_tiles_for_tms_with_non_standard_row_col_order():
         "+proj=s2 +lat_0=0.0 +lon_0=-90.0 +ellps=WGS84 +UVtoST=quadratic"
     )
     extent = [0.0, 0.0, 1.0, 1.0]
-    s2f4 = morecantile.TileMatrixSet.custom(extent, crs, id="S2F4")
+    with pytest.warns(NonWGS84GeographicCRS):
+        s2f4 = morecantile.TileMatrixSet.custom(extent, crs, id="S2F4")
     overlapping_tiles = s2f4.tiles(-100, 27, -95, 33, [6])
     assert len(list(overlapping_tiles)) == 30