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 342d621..7cf6b60 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,6 +1,34 @@ ## Unreleased +* 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": "http://www.opengis.net/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) diff --git a/morecantile/models.py b/morecantile/models.py index 9f5f9ab..aba911a 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): @@ -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,38 @@ 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 authority_code := feature_crs.to_authority(min_confidence=20): + authority, code = authority_code + feat.update( + { + "crs": { + "type": "name", + "properties": {"name": CRS_to_uri(feature_crs)}, + } } - } - ) + ) + else: + feat.update( + { + "crs": { + "type": "wkt", + "properties": {"wkt": feature_crs.to_wkt()}, + } + } + ) 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..aedf300 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 @@ -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' ) diff --git a/tests/test_morecantile.py b/tests/test_morecantile.py index 6530a72..0a059e5 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