From 64b4fb7e684f4b7c40ddd3e0ec38fb11e0cb7925 Mon Sep 17 00:00:00 2001 From: vincentsarago Date: Fri, 8 Nov 2024 15:02:10 +0100 Subject: [PATCH 1/5] use WGS84 as default CRS for GeoJSON outputs --- .gitignore | 1 + CHANGES.md | 5 +++++ morecantile/models.py | 41 +++++++++++++++++++++++++++++--------- morecantile/scripts/cli.py | 22 ++++++++++++++++++++ tests/test_cli.py | 29 +++++++++++++++++++++++++++ tests/test_morecantile.py | 25 +++++++++++++++++++++++ 6 files changed, 114 insertions(+), 9 deletions(-) diff --git a/.gitignore b/.gitignore index 788139f..fb17022 100644 --- a/.gitignore +++ b/.gitignore @@ -134,3 +134,4 @@ dmypy.json # VSCode .vscode .vscode/ +.notebooks/ diff --git a/CHANGES.md b/CHANGES.md index d82944c..9962dac 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,4 +1,9 @@ +## 6.2.0 (2024-11-08) + +* use `WGS84` as default CRS for GeoJSON output (as per specification) +* add `geographic_crs` option for `TileMatrixSet.feature` method + ## 6.1.0 (2024-10-17) * add `_tile_matrices_idx: Dict[str, int]` private attribute to improve `matrices` lookup diff --git a/morecantile/models.py b/morecantile/models.py index 84bcd08..6258858 100644 --- a/morecantile/models.py +++ b/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): @@ -1306,6 +1307,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. @@ -1327,16 +1329,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 ) @@ -1367,20 +1380,30 @@ def feature( }, } - 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, ) - 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) diff --git a/morecantile/scripts/cli.py b/morecantile/scripts/cli.py index ef6b309..c0070fb 100644 --- a/morecantile/scripts/cli.py +++ b/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 diff --git a/tests/test_cli.py b/tests/test_cli.py index 055b5d1..34ca6dc 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,5 +1,7 @@ """Tests of the morecantile CLI""" +import json + import pytest from click.testing import CliRunner @@ -20,6 +22,33 @@ def test_cli_shapes(): == '{"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' ) + 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": {"grid_crs": "http://www.opengis.net/def/crs/EPSG/0/3857", "grid_name": "WebMercatorQuad", "title": "XYZ tile (106, 193, 9)"}, "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 diff --git a/tests/test_morecantile.py b/tests/test_morecantile.py index c49f8cc..1e4b435 100644 --- a/tests/test_morecantile.py +++ b/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 From ad3c19cef60d1d468d499d98ba6494a130656cad Mon Sep 17 00:00:00 2001 From: vincentsarago Date: Thu, 19 Dec 2024 18:31:15 +0100 Subject: [PATCH 2/5] rename tms properties name --- morecantile/models.py | 4 ++-- tests/test_cli.py | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/morecantile/models.py b/morecantile/models.py index 5ce921e..c05ccd1 100644 --- a/morecantile/models.py +++ b/morecantile/models.py @@ -1377,8 +1377,8 @@ 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), }, } diff --git a/tests/test_cli.py b/tests/test_cli.py index 34ca6dc..aedf300 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -19,7 +19,7 @@ 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( @@ -28,7 +28,7 @@ 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' ) # With TMS's CRS @@ -54,7 +54,7 @@ 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' ) # buffer @@ -64,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 @@ -135,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' ) From 5526dd1dceac21ce46e84c94e7eed9bf287513e7 Mon Sep 17 00:00:00 2001 From: vincentsarago Date: Fri, 20 Dec 2024 14:10:14 +0100 Subject: [PATCH 3/5] use wkt --- morecantile/models.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/morecantile/models.py b/morecantile/models.py index c05ccd1..4321a38 100644 --- a/morecantile/models.py +++ b/morecantile/models.py @@ -1389,12 +1389,15 @@ def feature( UserWarning, stacklevel=1, ) - if epsg := feature_crs.to_epsg(): + if authority_code := feature_crs.to_authority(min_confidence=20): + authority, code = authority_code feat.update( { "crs": { - "type": "EPSG", - "properties": {"code": epsg}, + "type": "name", + "properties": { + "name": f"urn:ogc:def:crs:{authority}:0:{code}" + }, } } ) @@ -1402,8 +1405,8 @@ def feature( feat.update( { "crs": { - "type": "name", - "properties": {"name": CRS_to_uri(feature_crs)}, + "type": "wkt", + "properties": {"wkt": feature_crs.to_wkt()}, } } ) From e9c071d77c1006d476d46a354d789d9b8e9690f8 Mon Sep 17 00:00:00 2001 From: vincentsarago Date: Fri, 20 Dec 2024 14:18:47 +0100 Subject: [PATCH 4/5] update changelog --- CHANGES.md | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/CHANGES.md b/CHANGES.md index 985dc0f..b1c6d8b 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,8 +1,34 @@ ## Unreleased -* use `WGS84` as default CRS for GeoJSON output (as per specification) +* use `WGS84` as default CRS for `TileMatrixSet.feature` GeoJSON response (as per specification) * add `geographic_crs` option for `TileMatrixSet.feature` method +* update non-WGS84 CRS notation in `TileMatrixSet.feature` GeoJSON response + + ```python + # before + "properties": { + "crs": { + "type": "EPSG", + "properties": {"code": 3857}, + } + } + + # now + "properties": { + "crs": { + "type": "name", + "properties": {"name": f"urn:ogc:def:crs:EPSG:0:3857"}, + } + # or + "crs": { + "type": "wkt", + "properties": {"wkt": "..."}}, + } + } + ``` + +* rename `grid_name -> tms` and `grid_crs -> tms_crs` property names in `TileMatrixSet.feature` GeoJSON response * remove python 3.8 support ## 6.2.0 (2024-12-19) From d88897b88da390cc186977f8b58b882fa0fe2e7e Mon Sep 17 00:00:00 2001 From: vincentsarago Date: Fri, 20 Dec 2024 14:21:35 +0100 Subject: [PATCH 5/5] use CRS to URI --- CHANGES.md | 2 +- morecantile/models.py | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index b1c6d8b..7cf6b60 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -18,7 +18,7 @@ "properties": { "crs": { "type": "name", - "properties": {"name": f"urn:ogc:def:crs:EPSG:0:3857"}, + "properties": {"name": "http://www.opengis.net/def/crs/EPSG/0/3857"}, } # or "crs": { diff --git a/morecantile/models.py b/morecantile/models.py index 4321a38..aba911a 100644 --- a/morecantile/models.py +++ b/morecantile/models.py @@ -1389,15 +1389,14 @@ def feature( UserWarning, stacklevel=1, ) + if authority_code := feature_crs.to_authority(min_confidence=20): authority, code = authority_code feat.update( { "crs": { "type": "name", - "properties": { - "name": f"urn:ogc:def:crs:{authority}:0:{code}" - }, + "properties": {"name": CRS_to_uri(feature_crs)}, } } )