use WGS84 as default CRS for GeoJSON outputs

.gitignore

@@ -134,3 +134,4 @@ dmypy.json
 # VSCode
 .vscode
 .vscode/
+.notebooks/

CHANGES.md

@@ -1,6 +1,8 @@
 
 ## Unreleased
 
+* use `WGS84` as default CRS for GeoJSON output (as per specification)
+* add `geographic_crs` option for `TileMatrixSet.feature` method
 * remove python 3.8 support
 
 ## 6.2.0 (2024-12-19)

morecantile/models.py

@@ -41,6 +41,7 @@
 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)
 
 
 class CRSUri(BaseModel):
@@ -1308,6 +1309,7 @@ def feature(
         buffer: Optional[NumType] = None,
         precision: Optional[int] = None,
         projected: bool = False,
+        geographic_crs: Optional[CRS] = None,
     ) -> Dict:
         """
         Get the GeoJSON feature corresponding to a tile.
@@ -1329,16 +1331,27 @@ def feature(
             otherwise original coordinate values will be preserved (default).
         projected : bool, optional
             Return coordinates in TMS projection. Default is false.
+        geographic_crs: pyproj.CRS, optional
+            Geographic CRS to use when `projected=False`. Default to 'EPSG:4326' as per GeoJSON specification.
+            .
 
         Returns
         -------
         dict
 
         """
+        geographic_crs = geographic_crs or WGS84_CRS
+
+        feature_crs = self.crs._pyproj_crs
         west, south, east, north = self.xy_bounds(tile)
 
         if not projected:
-            west, south, east, north = self._to_geographic.transform_bounds(
+            feature_crs = geographic_crs
+            tr = pyproj.Transformer.from_crs(
+                self.crs._pyproj_crs, geographic_crs, always_xy=True
+            )
+
+            west, south, east, north = tr.transform_bounds(
                 west, south, east, north, densify_pts=21
             )
 
@@ -1364,26 +1377,36 @@ def feature(
             "geometry": geom,
             "properties": {
                 "title": f"XYZ tile {xyz}",
-                "grid_name": self.id,
-                "grid_crs": CRS_to_uri(self.crs._pyproj_crs),
+                "tms": self.id,
+                "tms_crs": CRS_to_uri(self.crs._pyproj_crs),
             },
         }
 
-        if projected:
+        if feature_crs != WGS84_CRS:
             warnings.warn(
                 "CRS is no longer part of the GeoJSON specification."
                 "Other projection than EPSG:4326 might not be supported.",
                 UserWarning,
                 stacklevel=1,
             )
-            feat.update(
-                {
-                    "crs": {
-                        "type": "EPSG",
-                        "properties": {"code": self.crs.to_epsg()},
+            if epsg := feature_crs.to_epsg():
+                feat.update(
+                    {
+                        "crs": {
+                            "type": "EPSG",
+                            "properties": {"code": epsg},
+                        }
                     }
-                }
-            )
+                )
+            else:
+                feat.update(
+                    {
+                        "crs": {
+                            "type": "name",
+                            "properties": {"name": CRS_to_uri(feature_crs)},
+                        }
+                    }
+                )
 
         if props:
             feat["properties"].update(props)

morecantile/scripts/cli.py

@@ -189,6 +189,11 @@ def cli(ctx, verbose, quiet):
     help="Path to TileMatrixSet JSON file.",
     type=click.Path(),
 )
+@click.option(
+    "--crs",
+    help="Geographic CRS. Default to WGS84.",
+    type=str,
+)
 @click.pass_context
 def shapes(  # noqa: C901
     ctx,
@@ -204,6 +209,7 @@ def shapes(  # noqa: C901
     extents,
     buffer,
     tms,
+    crs,
 ):
     """
     Reads one or more Web Mercator tile descriptions
@@ -223,6 +229,10 @@ def shapes(  # noqa: C901
     the properties object of the output feature.
 
     """
+    geographic_crs = None
+    if crs:
+        geographic_crs = CRS.from_user_input(crs)
+
     tilematrixset = morecantile.tms.get(identifier)
     if tms:
         with open(tms, "r") as f:
@@ -259,6 +269,7 @@ def shapes(  # noqa: C901
             projected=projected,
             buffer=buffer,
             precision=precision,
+            geographic_crs=geographic_crs,
         )
         bbox = feature["bbox"]
         w, s, e, n = bbox
@@ -526,6 +537,11 @@ def custom(
     default=None,
     help="Shift shape x and y values by a constant number",
 )
+@click.option(
+    "--crs",
+    help="Geographic CRS. Default to WGS84.",
+    type=str,
+)
 def tms_to_geojson(  # noqa: C901
     input,
     level,
@@ -538,8 +554,13 @@ def tms_to_geojson(  # noqa: C901
     collect,
     extents,
     buffer,
+    crs,
 ):
     """Print TMS document as GeoJSON."""
+    geographic_crs = None
+    if crs:
+        geographic_crs = CRS.from_user_input(crs)
+
     tms = morecantile.TileMatrixSet(**json.load(input))
     matrix = tms.matrix(level)
 
@@ -568,6 +589,7 @@ def tms_to_geojson(  # noqa: C901
                 projected=projected,
                 buffer=buffer,
                 precision=precision,
+                geographic_crs=geographic_crs,
             )
             bbox = feature["bbox"]
             w, s, e, n = bbox

tests/test_cli.py

@@ -1,5 +1,7 @@
 """Tests of the morecantile CLI"""
 
+import json
+
 import pytest
 from click.testing import CliRunner
 
@@ -17,15 +19,42 @@ def test_cli_shapes():
     assert result.exit_code == 0
     assert (
         result.output
-        == '{"bbox": [-105.46875, 39.909736, -104.765625, 40.446947], "geometry": {"coordinates": [[[-105.46875, 39.909736], [-105.46875, 40.446947], [-104.765625, 40.446947], [-104.765625, 39.909736], [-105.46875, 39.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"grid_crs": "http://www.opengis.net/def/crs/EPSG/0/3857", "grid_name": "WebMercatorQuad", "title": "XYZ tile (106, 193, 9)"}, "type": "Feature"}\n'
+        == '{"bbox": [-105.46875, 39.909736, -104.765625, 40.446947], "geometry": {"coordinates": [[[-105.46875, 39.909736], [-105.46875, 40.446947], [-104.765625, 40.446947], [-104.765625, 39.909736], [-105.46875, 39.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"title": "XYZ tile (106, 193, 9)", "tms": "WebMercatorQuad", "tms_crs": "http://www.opengis.net/def/crs/EPSG/0/3857"}, "type": "Feature"}\n'
     )
 
+    result = runner.invoke(
+        cli, ["shapes", "--precision", "6", "--geographic"], "[106, 193, 9]"
+    )
+    assert result.exit_code == 0
+    assert (
+        result.output
+        == '{"bbox": [-105.46875, 39.909736, -104.765625, 40.446947], "geometry": {"coordinates": [[[-105.46875, 39.909736], [-105.46875, 40.446947], [-104.765625, 40.446947], [-104.765625, 39.909736], [-105.46875, 39.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"title": "XYZ tile (106, 193, 9)", "tms": "WebMercatorQuad", "tms_crs": "http://www.opengis.net/def/crs/EPSG/0/3857"}, "type": "Feature"}\n'
+    )
+
+    # With TMS's CRS
+    with pytest.warns(UserWarning):
+        result = runner.invoke(
+            cli, ["shapes", "--precision", "6", "--projected"], "[106, 193, 9]"
+        )
+        assert result.exit_code == 0
+        feature = json.loads(result.output)
+        assert feature["crs"]
+
+    # geographic CRS (non WGS84)
+    with pytest.warns(UserWarning):
+        result = runner.invoke(
+            cli, ["shapes", "--precision", "6", "--crs", "epsg:4150"], "[106, 193, 9]"
+        )
+        assert result.exit_code == 0
+        feature = json.loads(result.output)
+        assert feature["crs"]
+
     # tile as arg
     result = runner.invoke(cli, ["shapes", "[106, 193, 9]", "--precision", "6"])
     assert result.exit_code == 0
     assert (
         result.output
-        == '{"bbox": [-105.46875, 39.909736, -104.765625, 40.446947], "geometry": {"coordinates": [[[-105.46875, 39.909736], [-105.46875, 40.446947], [-104.765625, 40.446947], [-104.765625, 39.909736], [-105.46875, 39.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"grid_crs": "http://www.opengis.net/def/crs/EPSG/0/3857", "grid_name": "WebMercatorQuad", "title": "XYZ tile (106, 193, 9)"}, "type": "Feature"}\n'
+        == '{"bbox": [-105.46875, 39.909736, -104.765625, 40.446947], "geometry": {"coordinates": [[[-105.46875, 39.909736], [-105.46875, 40.446947], [-104.765625, 40.446947], [-104.765625, 39.909736], [-105.46875, 39.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"title": "XYZ tile (106, 193, 9)", "tms": "WebMercatorQuad", "tms_crs": "http://www.opengis.net/def/crs/EPSG/0/3857"}, "type": "Feature"}\n'
     )
 
     # buffer
@@ -35,7 +64,7 @@ def test_cli_shapes():
     assert result.exit_code == 0
     assert (
         result.output
-        == '{"bbox": [-106.46875, 38.909736, -103.765625, 41.446947], "geometry": {"coordinates": [[[-106.46875, 38.909736], [-106.46875, 41.446947], [-103.765625, 41.446947], [-103.765625, 38.909736], [-106.46875, 38.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"grid_crs": "http://www.opengis.net/def/crs/EPSG/0/3857", "grid_name": "WebMercatorQuad", "title": "XYZ tile (106, 193, 9)"}, "type": "Feature"}\n'
+        == '{"bbox": [-106.46875, 38.909736, -103.765625, 41.446947], "geometry": {"coordinates": [[[-106.46875, 38.909736], [-106.46875, 41.446947], [-103.765625, 41.446947], [-103.765625, 38.909736], [-106.46875, 38.909736]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"title": "XYZ tile (106, 193, 9)", "tms": "WebMercatorQuad", "tms_crs": "http://www.opengis.net/def/crs/EPSG/0/3857"}, "type": "Feature"}\n'
     )
 
     # Output is compact
@@ -106,7 +135,7 @@ def test_cli_shapesWGS84():
     assert result.exit_code == 0
     assert (
         result.output
-        == '{"bbox": [-105.46875, 40.099155, -104.765625, 40.636956], "geometry": {"coordinates": [[[-105.46875, 40.099155], [-105.46875, 40.636956], [-104.765625, 40.636956], [-104.765625, 40.099155], [-105.46875, 40.099155]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"grid_crs": "http://www.opengis.net/def/crs/EPSG/0/3395", "grid_name": "WorldMercatorWGS84Quad", "title": "XYZ tile (106, 193, 9)"}, "type": "Feature"}\n'
+        == '{"bbox": [-105.46875, 40.099155, -104.765625, 40.636956], "geometry": {"coordinates": [[[-105.46875, 40.099155], [-105.46875, 40.636956], [-104.765625, 40.636956], [-104.765625, 40.099155], [-105.46875, 40.099155]]], "type": "Polygon"}, "id": "(106, 193, 9)", "properties": {"title": "XYZ tile (106, 193, 9)", "tms": "WorldMercatorWGS84Quad", "tms_crs": "http://www.opengis.net/def/crs/EPSG/0/3395"}, "type": "Feature"}\n'
     )
 
 

tests/test_morecantile.py

@@ -183,11 +183,36 @@ def test_feature():
         feat = tms.feature(
             morecantile.Tile(1, 0, 1), projected=True, fid="1", props={"some": "thing"}
         )
+    assert feat["crs"]
     assert feat["bbox"]
     assert feat["id"] == "1"
     assert feat["geometry"]
     assert len(feat["properties"].keys()) == 4
 
+    # These extent coordinates are in EPSG:2056 (CH)
+    custom_tms = morecantile.TileMatrixSet.custom(
+        [2696082.04374708, 1289407.53195196, 2696210.04374708, 1289535.53195196],
+        CRS.from_epsg("2056"),
+    )
+    assert custom_tms.geographic_crs != CRS.from_epsg(4326)
+    # Warn when geographic CRS is not WGS84
+    with pytest.warns(UserWarning):
+        feat = custom_tms.feature(
+            morecantile.Tile(1, 0, 1),
+            projected=False,
+            geographic_crs=custom_tms.geographic_crs,
+        )
+        assert feat["crs"]
+
+    # By default we use WGS84 CRS (as per GeoJSON spec)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        feat = custom_tms.feature(
+            morecantile.Tile(1, 0, 1),
+            projected=False,
+        )
+        assert not feat.get("crs")
+
 
 ################################################################################
 # replicate mercantile tests