diff --git a/.gitignore b/.gitignore index 788139f..223a84f 100644 --- a/.gitignore +++ b/.gitignore @@ -134,3 +134,5 @@ dmypy.json # VSCode .vscode .vscode/ + +.notebooks/ diff --git a/morecantile/errors.py b/morecantile/errors.py index 83af4be..73e0339 100644 --- 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 index 84bcd08..842ebf5 100644 --- 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 index 63d7bb0..fb944e0 100644 --- 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 index c1a84ca..d843423 100644 --- 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 index c49f8cc..a14616d 100644 --- 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