Skip to content

Commit

Permalink
posible fix for always assuming WGS84
Browse files Browse the repository at this point in the history
  • Loading branch information
AndrewAnnex committed Sep 30, 2024
1 parent 1829fe1 commit 02cefc8
Show file tree
Hide file tree
Showing 3 changed files with 32 additions and 9 deletions.
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -130,3 +130,6 @@ dmypy.json

# PyCharm:
.idea

# vscode:
.vscode
16 changes: 10 additions & 6 deletions morecantile/models.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,6 @@
NumType = Union[float, int]
BoundsType = Tuple[NumType, NumType]
LL_EPSILON = 1e-11
WGS84_CRS = pyproj.CRS.from_epsg(4326)
axesInfo = Annotated[List[str], Field(min_length=2, max_length=2)]


Expand Down Expand Up @@ -486,16 +485,18 @@ class TileMatrixSet(BaseModel, arbitrary_types_allowed=True):
]

# Private attributes
_geographic_crs: pyproj.CRS = PrivateAttr(default=WGS84_CRS)
_geographic_crs: pyproj.CRS = PrivateAttr()
_to_geographic: pyproj.Transformer = PrivateAttr()
_from_geographic: pyproj.Transformer = PrivateAttr()

def __init__(self, **data):
"""Set private attributes."""
super().__init__(**data)

# TODO is the geodetic_crs always the geographic crs?
# maybe instead switch depending on if it's a projection or a geodcrs/geogcrs
self._geographic_crs = pyproj.CRS.from_user_input(
data.get("_geographic_crs", WGS84_CRS)
data.get("_geographic_crs", self.crs._pyproj_crs.geodetic_crs)
)

try:
Expand Down Expand Up @@ -656,7 +657,7 @@ def custom(
title: Optional[str] = None,
id: Optional[str] = None,
ordered_axes: Optional[List[str]] = None,
geographic_crs: pyproj.CRS = WGS84_CRS,
geographic_crs: pyproj.CRS = None,
screen_pixel_size: float = 0.28e-3,
decimation_base: int = 2,
**kwargs: Any,
Expand Down Expand Up @@ -689,8 +690,8 @@ def custom(
Tile Matrix Set title
id: str, optional
Tile Matrix Set identifier
geographic_crs: pyproj.CRS
Geographic (lat,lon) coordinate reference system (default is EPSG:4326)
geographic_crs: pyproj.CRS, optional
Geographic (lat,lon) coordinate reference system (default is CRS's geodetic CRS)
ordered_axes: list of str, optional
Override Axis order (e.g `["N", "S"]`) else default to CRS's metadata
screen_pixel_size: float, optional
Expand All @@ -705,6 +706,9 @@ def custom(
TileMatrixSet
"""
# TODO is this correct to always assume?
geographic_crs = geographic_crs or crs.geodetic_crs

matrix_scale = matrix_scale or [1, 1]

if ordered_axes:
Expand Down
22 changes: 19 additions & 3 deletions tests/test_models.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from collections.abc import Iterable

import pyproj
import pyproj.crs
import pytest
from pydantic import ValidationError
from rasterio.crs import CRS as rioCRS
Expand Down Expand Up @@ -288,6 +289,19 @@ def test_zoom_for_res():
assert tms.zoom_for_res(10) == 14
assert tms.zoom_for_res(5000) == 6

def test_custom_not_earth():
crs = pyproj.CRS.from_user_input('IAU_2015:49900')
extent = [-90, -180, 90, 180]
mars_tms = morecantile.TileMatrixSet.custom(extent, crs, id="MarsGeographicCRS")
assert '4326' not in mars_tms.geographic_crs.to_wkt()
assert '4326' not in mars_tms.rasterio_geographic_crs.to_wkt()
wkt_mars_web_mercator = "PROJCRS[\"Mars (2015) - Sphere XY / Pseudo-Mercator\",BASEGEOGCRS[\"Mars (2015) - Sphere\",DATUM[\"Mars (2015) - Sphere\",ELLIPSOID[\"Mars (2015) - Sphere\",3396190,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],ANCHOR[\"Viking 1 lander : 47.95137 W\"]],PRIMEM[\"Reference Meridian\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CONVERSION[\"Popular Visualisation Pseudo-Mercator\",METHOD[\"Popular Visualisation Pseudo Mercator\",ID[\"EPSG\",1024]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"easting (X)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"northing (Y)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],USAGE[SCOPE[\"Web mapping and visualisation.\"],AREA[\"World between 85.06 S and 85.06 N.\"],BBOX[-85.850511287,-180,85.0511287,180]],REMARK[\"Use semi-major radius as sphere radius for interoperability. Source of IAU Coordinate systems: doi:10.1007/s10569-017-9805-5\"]]"
crs_mars_web_mercator = pyproj.CRS.from_wkt(wkt_mars_web_mercator)
extent_wm = [-85.850511287,-180,85.0511287,180]
mars_tms_wm = morecantile.TileMatrixSet.custom(extent_wm, 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()


def test_schema():
"""Translate Model to Schema."""
Expand All @@ -300,13 +314,15 @@ 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]
with pytest.warns(UserWarning):
tms = morecantile.TileMatrixSet.custom(extent, crs, id="MarsNPolek2MOLA5k")
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))
assert json_doc["crs"] == "http://www.opengis.net/def/crs/IAU/2015/49930"

# test to ensure a correct geographic crs
mars_tms = TileMatrixSet.model_validate(json_doc)
assert '4326' not in mars_tms.geographic_crs.to_wkt()
assert '4326' not in mars_tms.rasterio_geographic_crs.to_wkt()
crs = pyproj.CRS.from_epsg(3031)
extent = [-948.75, -543592.47, 5817.41, -3333128.95] # From https:///epsg.io/3031
tms = morecantile.TileMatrixSet.custom(extent, crs, id="MyCustomTmsEPSG3031")
Expand Down

0 comments on commit 02cefc8

Please sign in to comment.