diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..6062365 --- /dev/null +++ b/.flake8 @@ -0,0 +1,5 @@ +[flake8] +max-line-length = 88 +ignore = E203, E266, E501, W503, F403, F401, E741 +max-complexity = 18 +select = B,C,E,F,W,T4,B9 \ No newline at end of file diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml new file mode 100644 index 0000000..2705736 --- /dev/null +++ b/.github/workflows/python-package.yml @@ -0,0 +1,54 @@ +# This workflow will install Python dependencies, run tests and lint with a variety of Python versions +# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions + +name: Python package test + +on: + push: + branches: [main] + pull_request: + branches: [main] + +jobs: + test: + runs-on: ${{ matrix.os }} + + strategy: + fail-fast: false + matrix: + python-version: ["3.7", "3.8", "3.9", "3.10"] + os: ["ubuntu-22.04"] + + steps: + - uses: actions/checkout@v2 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + env: + CURL_CA_BUNDLE: /etc/ssl/certs/ca-certificates.crt + run: | + python -m pip install --upgrade pip + pip install -e .[test] + + - name: Lint with flake8 + run: | + # stop the build if there are Python syntax errors or undefined names + flake8 tilematrix --count --select=E9,F63,F7,F82 --show-source --statistics + # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide + flake8 tilematrix --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + + - name: Test with pytest + env: + CURL_CA_BUNDLE: /etc/ssl/certs/ca-certificates.crt + run: | + pytest -v --cov tilematrix + + - name: Coveralls + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + coveralls --service=github \ No newline at end of file diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml new file mode 100644 index 0000000..6d16234 --- /dev/null +++ b/.github/workflows/python-publish.yml @@ -0,0 +1,31 @@ +# This workflows will upload a Python Package using Twine when a release is created +# For more information see: https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#publishing-to-package-registries + +name: Upload Python Package + +on: + release: + types: [created] + +jobs: + deploy: + + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v2 + - name: Set up Python + uses: actions/setup-python@v2 + with: + python-version: '3.x' + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install hatch + - name: Build and publish + env: + HATCH_INDEX_USER: ${{ secrets.PYPI_USERNAME }} + HATCH_INDEX_AUTH: ${{ secrets.PYPI_PASSWORD }} + run: | + hatch build + hatch publish diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..198d4cf --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +*.pyc +dist/ +build/ +testenv/ +*.egg-info +test/testdata/out/* +venv/ +.cache +.eggs +.coverage +htmlcov +.pytest* \ No newline at end of file diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..8388ca1 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,24 @@ +repos: +- repo: https://github.com/ambv/black + rev: 23.3.0 + hooks: + - id: black + language_version: python3 +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v1.2.3 + hooks: + - id: flake8 +- repo: https://github.com/PyCQA/flake8 + rev: 6.0.0 + hooks: + - id: flake8 +- repo: https://github.com/PyCQA/autoflake + rev: v2.1.1 + hooks: + - id: autoflake +- repo: https://github.com/pycqa/isort + rev: 5.12.0 + hooks: + - id: isort + name: isort (python) + args: ["--profile", "black"] diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..306f58e --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,16 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "Python: Current File", + "type": "python", + "request": "launch", + "program": "${file}", + "console": "integratedTerminal", + "justMyCode": true + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..b2b8866 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,7 @@ +{ + "python.testing.pytestArgs": [ + "test" + ], + "python.testing.unittestEnabled": false, + "python.testing.pytestEnabled": true +} \ No newline at end of file diff --git a/CHANGELOG.rst b/CHANGELOG.rst new file mode 100644 index 0000000..bf099ae --- /dev/null +++ b/CHANGELOG.rst @@ -0,0 +1,214 @@ +######### +Changelog +######### + +2023.12.0 - 2022-12-15 +---------------------- +* fix grid warning +* add more pre-commit hooks + + +2022.12.0 - 2022-12-16 +---------------------- +* make fit for `Shapely 2.0` + + +2022.11.3 - 2022-11-08 +---------------------- +* replace `setuptools` with `hatch` + + +2022.11.2 - 2022-11-08 +---------------------- +* remove `conda` build files +* `tilematrix-feedstock `_ `conda-forge feedstock` repository for releasing new versions on `conda-forge` + + +2022.11.1 - 2022-11-03 +---------------------- +* fix `conda` & `pip ` builds + + +2022.11.0 - 2022-11-03 +---------------------- +* test also for python `3.9` +* renaming `master` branch to `main` +* add ``conda`` publish to github actions (workflows) + + +2022.3.0 - 2022-03-15 +--------------------- +* add option to exactly get intersection tiles using `TilePyramid.tiles_from_geom(exact=True)` + + +2021.11.0 - 2021-11-12 +---------------------- +* enable yielding tiles in batches by either row or column for following methods: + * `TilePyramid.tiles_from_bounds()` + * `TilePyramid.tiles_from_bbox()` + * `TilePyramid.tiles_from_geom()` + +* convert TilePyramid arguments into keyword arguments + + +0.21 +---- +* allow metatiling up to 512 +* use GitHub actions instead of travis +* use black and flake8 pre-commit checks + + +0.20 +---- +* fixed pixel size calculation on irregular grids with metatiling (closing #33) +* ``TilePyramid.tile_x_size()``, ``TilePyramid.tile_y_size()``, ``TilePyramid.tile_height()``, ``TilePyramid.tile_width()`` are deprecated +* metatiles are clipped to ``TilePyramid.bounds`` but ``pixelbuffer`` of edge tiles can exceed them unless it is a global grid + +0.19 +---- +* Python 2 not supported anymore +* ``TilePyramid.srid`` and ``TilePyramid.type`` are deprecated +* ``GridDefinition`` can now be loaded from package root +* ``GridDefinition`` got ``to_dict()`` and ``from_dict()`` methods + + +0.18 +---- +* order of ``Tile.shape`` swapped to ``(height, width)`` in order to match rasterio array interpretation + +0.17 +---- +* make ``Tile`` iterable to enable ``tuple(Tile)`` return the tile index as tuple + +0.16 +---- +* make ``Tile`` objects hashable & comparable + +0.15 +---- +* add ``snap_bounds()`` function +* add ``snap_bounds`` command to ``tmx`` +* add ``snap_bbox`` command to ``tmx`` +* in ``tile_from_xy()`` add ``on_edge_use`` option specify behavior when point hits grid edges +* cleaned up ``_tiles_from_cleaned_bounds()`` and ``tile_from_xy()`` functions + +------ +0.14.2 +------ +* attempt to fix ``tmx`` command when installing tilematrix via pip + +---- +0.14 +---- +* add ``tmx`` CLI with subcommands: + * `bounds`: Print bounds of given Tile. + * `bbox`: Print bounding box geometry of given Tile. + * `tile`: Print Tile covering given point. + * `tiles`: Print Tiles covering given bounds. + +---- +0.13 +---- +* fixed ``tiles_from_geom()`` bug when passing on a Point (#19) +* add ``tile_from_xy()`` function + +---- +0.12 +---- +* added better string representations for ``Tile`` and ``TilePyramid`` +* added ``GridDefinition`` to better handle custom grid parameters +* ``TilePyramid`` instances are now comparable by ``==`` and ``!=`` + +---- +0.11 +---- +* custom grid defnitions enabled + +--- +0.10 +--- +* new tag for last version to fix Python 3 build + +--- +0.9 +--- +* added Python 3 support +* use NamedTuple for Tile index + +--- +0.8 +--- +* ``intersecting`` function fixed (rounding error caused return of wrong tiles) + +--- +0.7 +--- +* converted tuples for bounds and shape attributes to namedtuples + +--- +0.6 +--- +* added ``pytest`` and test cases +* fixed metatiling shape error on low zoom levels +* split up code into internal modules +* travis CI and coveralls.io integration + +--- +0.5 +--- +* ``intersection()`` doesn't return invalid tiles. +* Moved copyright to EOX IT Services + +--- +0.4 +--- +* Decision to remove ``MetaTilePyramid`` class (now returns a ``DeprecationWarning``). +* TilePyramid now has its own ``metatiling`` parameter. +* ``intersecting()`` function for ``Tile`` and ``TilePyramid`` to relate between ``TilePyramids`` with different ``metatiling`` settings. + +--- +0.3 +--- +* fixed duplicate tile return in tiles_from_bounds() +* rasterio's CRS() class replaced CRS dict + +--- +0.2 +--- +* introduced handling of antimeridian: + * ``get_neighbor()`` also gets tiles from other side + * ``.shape()`` returns clipped tile shape + * added ``tiles_from_bounds()`` + * added ``clip_geometry_to_srs_bounds()`` + +--- +0.1 +--- +* added Spherical Mercator support +* removed IO module (moved to `mapchete `_) +* removed deprecated ``OutputFormats`` +* introduced ``get_parent()`` and ``get_children()`` functions for ``Tile`` + +----- +0.0.4 +----- +* introduced ``Tile`` object +* read_raster_window() is now a generator which returns only a numpy array +* read_vector_window() is a generator which returns a GeoJSON-like object with a geometry clipped to tile boundaries +* proper error handling (removed ``sys.exit(0)``) + +----- +0.0.3 +----- +* rewrote io module +* separated and enhanced OutputFormats + +----- +0.0.2 +----- +* fixed wrong link to github repository + +----- +0.0.1 +----- +* basic functionality diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..f88dd0f --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2015 - 2019 EOX IT Services GmbH + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.rst b/README.rst new file mode 100644 index 0000000..8b425af --- /dev/null +++ b/README.rst @@ -0,0 +1,142 @@ +========== +TilePyramidCut +========== + +TilePyramidCut removes files from a tile pyramids. + +.. image:: https://badge.fury.io/py/tilematrix.svg + :target: https://badge.fury.io/py/tilematrix + +.. image:: https://travis-ci.org/ungarj/tilematrix.svg?branch=master + :target: https://travis-ci.org/ungarj/tilematrix + +.. image:: https://coveralls.io/repos/github/ungarj/tilematrix/badge.svg?branch=master + :target: https://coveralls.io/github/ungarj/tilematrix?branch=master + +.. image:: https://anaconda.org/conda-forge/tilematrix/badges/version.svg + :target: https://anaconda.org/conda-forge/tilematrix + +.. image:: https://img.shields.io/pypi/pyversions/mapchete.svg + + +The module is designed to translate between tile indices (zoom, row, column = ZYX) and +map coordinates (e.g. latitute, longitude). + +Tilematrix supports **metatiling** and **tile buffers**. Furthermore it makes +heavy use of shapely_ and it can also generate ``Affine`` objects per tile which +facilitates working with rasterio_ for tile based data reading and writing. + +It is very similar to mercantile_ but besides of supporting spherical mercator +tile pyramids, it also supports geodetic (WGS84) tile pyramids. + +.. _shapely: http://toblerity.org/shapely/ +.. _rasterio: https://github.com/mapbox/rasterio +.. _mercantile: https://github.com/mapbox/mercantile + +------------ +Installation +------------ + +Use ``conda`` to install the latest stable version: + +.. code-block:: shell + + conda install -c conda-forge -y tilematrix + +Use ``pip`` to install the latest stable version: + +.. code-block:: shell + + pip install tilematrix + +Manually install the latest development version + +.. code-block:: shell + + pip install -r requirements.txt + python setup.py install + + +------------- +Documentation +------------- + +* `API documentation `_ +* `examples `_ + +CLI +--- + +This package ships with a command line tool ``tmx`` which provides the following +subcommands: + +* ``bounds``: Print bounds of given Tile. +* ``bbox``: Print bounding box geometry of given Tile. +* ``tile``: Print Tile covering given point. +* ``tiles``: Print Tiles covering given bounds. + +Geometry outputs can either be formatted as ``WKT`` or ``GeoJSON``. For example +the following command will print a valid ``GeoJSON`` representing all tiles +for zoom level 1 of the ``geodetic`` WMTS grid: + +.. code-block:: shell + + $ tmx -f GeoJSON tiles -- 1 -180 -90 180 90 + { + "type": "FeatureCollection", + "features": [ + {"geometry": {"coordinates": [[[-90.0, 0.0], [-90.0, 90.0], [-180.0, 90.0], [-180.0, 0.0], [-90.0, 0.0]]], "type": "Polygon"}, "properties": {"col": 0, "row": 0, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[0.0, 0.0], [0.0, 90.0], [-90.0, 90.0], [-90.0, 0.0], [0.0, 0.0]]], "type": "Polygon"}, "properties": {"col": 1, "row": 0, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[90.0, 0.0], [90.0, 90.0], [0.0, 90.0], [0.0, 0.0], [90.0, 0.0]]], "type": "Polygon"}, "properties": {"col": 2, "row": 0, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[180.0, 0.0], [180.0, 90.0], [90.0, 90.0], [90.0, 0.0], [180.0, 0.0]]], "type": "Polygon"}, "properties": {"col": 3, "row": 0, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[-90.0, -90.0], [-90.0, 0.0], [-180.0, 0.0], [-180.0, -90.0], [-90.0, -90.0]]], "type": "Polygon"}, "properties": {"col": 0, "row": 1, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[0.0, -90.0], [0.0, 0.0], [-90.0, 0.0], [-90.0, -90.0], [0.0, -90.0]]], "type": "Polygon"}, "properties": {"col": 1, "row": 1, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[90.0, -90.0], [90.0, 0.0], [0.0, 0.0], [0.0, -90.0], [90.0, -90.0]]], "type": "Polygon"}, "properties": {"col": 2, "row": 1, "zoom": 1}, "type": "Feature"}, + {"geometry": {"coordinates": [[[180.0, -90.0], [180.0, 0.0], [90.0, 0.0], [90.0, -90.0], [180.0, -90.0]]], "type": "Polygon"}, "properties": {"col": 3, "row": 1, "zoom": 1}, "type": "Feature"} + ] + } + + + +Print ``WKT`` representation of tile ``4 15 23``: + +.. code-block:: shell + + $ tmx bbox 4 15 23 + POLYGON ((90 -90, 90 -78.75, 78.75 -78.75, 78.75 -90, 90 -90)) + + +Also, tiles can have buffers around called ``pixelbuffer``: + +.. code-block:: shell + + $ tmx --pixelbuffer 10 bbox 4 15 23 + POLYGON ((90.439453125 -90, 90.439453125 -78.310546875, 78.310546875 -78.310546875, 78.310546875 -90, 90.439453125 -90)) + + +Print ``GeoJSON`` representation of tile ``4 15 23`` on a ``mercator`` tile +pyramid: + +.. code-block:: shell + + $ tmx -output_format GeoJSON -grid mercator bbox 4 15 15 + {"type": "Polygon", "coordinates": [[[20037508.342789203, -20037508.3427892], [20037508.342789203, -17532819.799940553], [17532819.799940553, -17532819.799940553], [17532819.799940553, -20037508.3427892], [20037508.342789203, -20037508.3427892]]]} + +---------------- +Conda Publishing +---------------- + +Use bot pull requests generated with every release at tilematrix-feedstock_ repository for releasing new versions on ``conda-forge``. + + +------- +License +------- + +MIT License + +Copyright (c) 2015-2022 `EOX IT Services`_ + +.. _`EOX IT Services`: https://eox.at/ + +.. _`tilematrix-feedstock`: https://github.com/conda-forge/tilematrix-feedstock diff --git a/doc/examples.md b/doc/examples.md new file mode 100644 index 0000000..922ffb6 --- /dev/null +++ b/doc/examples.md @@ -0,0 +1,85 @@ +# Examples + +```python +from tilematrix import TilePyramid +from shapely.geometry import Polygon + +# initialize TilePyramid +tile_pyramid = TilePyramid("geodetic") + +some_polygon = Polygon([(0, 0), (1, 1), (1, 0)]) +zoom = 8 + +``` + +Now, let's get all tile IDs (zoom, row, col) of tiles intersecting with our +example geometry: + +```python +for tile in tile_pyramid.tiles_from_geom(some_polygon, zoom): + print tile.id + +``` + +output: +```python +(8, 126, 256) +(8, 126, 257) +(8, 127, 256) +(8, 127, 257) + +``` + +Use the ``bbox()`` method to get the tile bounding box geometries: +```python +for tile in tile_pyramid.tiles_from_geom(some_polygon, zoom): + print tile.bbox() +``` + +output: +```python +POLYGON ((0 1.40625, 0.703125 1.40625, 0.703125 0.703125, 0 0.703125, 0 1.40625)) +POLYGON ((0.703125 1.40625, 1.40625 1.40625, 1.40625 0.703125, 0.703125 0.703125, 0.703125 1.40625)) +POLYGON ((0 0.703125, 0.703125 0.703125, 0.703125 0, 0 0, 0 0.703125)) +POLYGON ((0.703125 0.703125, 1.40625 0.703125, 1.40625 0, 0.703125 0, 0.703125 0.703125)) +``` + +We can also create a buffer around the tiles, using the ``pixelbuffer`` argument. A value of 1 will create a 1 px buffer around the tile. This depends on the initial pixel "resolution" the ``TilePyramid`` was initialized with (default ``256``). +```python +for tile in tile_pyramid.tiles_from_geom(some_polygon, zoom): + print tile.bbox(pixelbuffer=1) +``` + +output: +```python +POLYGON ((-0.00274658203125 1.40899658203125, 0.70587158203125 1.40899658203125, 0.70587158203125 0.70037841796875, -0.00274658203125 0.70037841796875, -0.00274658203125 1.40899658203125)) +POLYGON ((0.70037841796875 1.40899658203125, 1.40899658203125 1.40899658203125, 1.40899658203125 0.70037841796875, 0.70037841796875 0.70037841796875, 0.70037841796875 1.40899658203125)) +POLYGON ((-0.00274658203125 0.70587158203125, 0.70587158203125 0.70587158203125, 0.70587158203125 -0.00274658203125, -0.00274658203125 -0.00274658203125, -0.00274658203125 0.70587158203125)) +POLYGON ((0.70037841796875 0.70587158203125, 1.40899658203125 0.70587158203125, 1.40899658203125 -0.00274658203125, 0.70037841796875 -0.00274658203125, 0.70037841796875 0.70587158203125)) +``` + +## Metatiling + +Metatiling is often used to combine smaller ``(256, 256)`` tiles into bigger ones. The metatiling parameter describes how many small tiles are combines into one bigger metatile. For example a metatiling parameter of 2 would combine 2x2 smaller tiles into one metatile. + +You can activate metatiling by initializing a ``TilePyramid`` with a metatiling value. + +```python +# initialize TilePyramid +tile_pyramid = TilePyramid("geodetic", metatiling=2) +``` + +As the tiles are now bigger, the code we used above: +```python +some_polygon = Polygon([(0, 0), (1, 1), (1, 0)]) +zoom = 8 + +for tile in tile_pyramid.tiles_from_geom(some_polygon, zoom): + print tile.bbox() +``` + +now returns just one but bigger metatile: +output: +```python +POLYGON ((1.40625 0, 1.40625 1.40625, 0 1.40625, 0 0, 1.40625 0)) +``` diff --git a/doc/tilematrix.md b/doc/tilematrix.md new file mode 100644 index 0000000..c6ee5e1 --- /dev/null +++ b/doc/tilematrix.md @@ -0,0 +1,151 @@ +# tilematrix + +The two base classes are ``TilePyramid`` which defines tile matrices in various zoom levels and its members, the``Tile`` objects. ``TilePyramid`` + +## TilePyramid + +```python +TilePyramid(grid_definition, tile_size=256) +``` +* ``grid_definition``: Either one of the predefined grids (``geodetic`` or + ``mercator``) or a custom grid definition in form of a dictionary. For example: + ```python + { + "shape": (1, 1), + "bounds": (2426378.0132, 1528101.2618, 6293974.6215, 5395697.8701), + "is_global": False, + "epsg": 3035 + } + ``` + * ``shape`` (height, width): Indicates the number of tiles per column and row at **zoom level 0**. + * ``bounds`` (left, bottom, right, top): Units are CRS units. + + Please note that the aspect ratios of ``shape`` and ``bounds`` have to be the same. + * Alternatively to ``epsg``, a custom ``proj`` string can be used: + ```python + "proj": """ + +proj=ortho + +lat_0=-90 + +lon_0=0 + +x_0=0 + +y_0=0 + +ellps=WGS84 + +units=m +no_defs + """ + ``` + * ``is_global``: Indicates whether the grid covers the whole globe or just a region. + +* ``tile_size``: Optional, specifies the target resolution of each tile (i.e. each tile will have 256x256 px). Default is ``256``. + +### Variables +* ``type``: The projection it was initialized with (``geodetic``, ``mercator`` or ``custom``). +* ``tile_size``: The pixelsize per tile it was initialized with (default ``256``). +* ``left``: Left boundary of tile matrix. +* ``top``: Top boundary of tile matrix. +* ``right ``: Right boundary of tile matrix. +* ``bottom ``: Bottom boundary of tile matrix. +* ``x_size``: Horizontal size of tile matrix in map coordinates. +* ``y_size``: Vertical size of tile matrix in map coordinates. +* ``crs``: ``rasterio.crs.CRS()`` object. +* ``srid``: Spatial reference ID (EPSG code) if available, else ``None``. + +### Methods + +#### Matrix properties +* ``matrix_width(zoom)``: Returns the number of columns in the current tile matrix. + * ``zoom``: Zoom level. +* ``matrix_height(zoom)``: Returns the number of rows in the current tile matrix. + * ``zoom``: Zoom level. + +#### Geometry +* ``tile_x_size(zoom)``: Returns a float, indicating the width of each tile at this zoom level in ``TilePyramid`` CRS units. + * ``zoom``: Zoom level. +* ``tile_y_size(zoom)``: Returns a floats, indicating the height of each tile at this zoom level in ``TilePyramid`` CRS units. + * ``zoom``: Zoom level. +* ``pixel_x_size(zoom)``: Returns a float, indicating the vertical pixel size in ``TilePyramid`` CRS units. + * ``zoom``: Zoom level. +* ``pixel_y_size(zoom)``: Returns a float, indicating the horizontal pixel size in ``TilePyramid`` CRS units. + * ``zoom``: Zoom level. +* ``top_left_tile_coords(zoom, row, col)``: Returns a tuple of two floats, indicating the top left coordinates of the given tile. + * ``zoom``: Zoom level. + * ``row``: Row in ``TilePyramid``. + * ``col``: Column in ``TilePyramid``. +* ``tile_bounds(zoom, row, col, pixelbuffer=None)``: Returns ``left``, ``bottom``, ``right``, ``top`` coordinates of given tile. + * ``zoom``: Zoom level. + * ``row``: Row in ``TilePyramid``. + * ``col``: Column in ``TilePyramid``. + * ``pixelbuffer``: Optional buffer around tile boundaries in pixels. +* ``tile_bbox(zoom, row, col, pixelbuffer=None)``: Returns the bounding box for given tile as a ``Polygon``. + * ``zoom``: Zoom level. + * ``row``: Row in ``TilePyramid``. + * ``col``: Column in ``TilePyramid``. + * ``pixelbuffer``: Optional buffer around tile boundaries in pixels. + + +#### Tiles +* ``intersecting(tile)``: Return all tiles intersecting with tile. This helps translating between TilePyramids with different metatiling settings. + * ``tile``: ``Tile`` object. +* ``tiles_from_bbox(geometry, zoom)``: Returns tiles intersecting with the given bounding box at given zoom level. + * ``geometry``: Must be a ``Polygon`` object. + * ``zoom``: Zoom level. +* ``tiles_from_geom(geometry, zoom)``: Returns tiles intersecting with the given geometry at given zoom level. + * ``geometry``: Must be one out of ``Polygon``, ``MultiPolygon``, ``LineString``, ``MultiLineString``, ``Point``, ``MultiPoint``. + * ``zoom``: Zoom level. + + +## Tile + +```python +Tile(tile_pyramid, zoom, row, col) +``` +* ``tile_pyramid``: A ``TilePyramid`` this Tile belongs to. +* ``zoom``: Tile zoom level. +* ``row``: Tile row. +* ``col``: Tile column. + +### Variables +After initializing, the object has the following properties: +* ``tile_pyramid`` or ``tp``: ``TilePyramid`` this Tile belongs to. +* ``crs``: ``rasterio.crs.CRS()`` object. +* ``srid``: Spatial reference ID (EPSG code) if available, else ``None``. +* ``zoom``: Zoom. +* ``row``: Row. +* ``col``: Col. +* ``x_size``: Horizontal size in SRID units. +* ``y_size``: Vertical size in SRID units. +* ``index`` or ``id``: Tuple of ``(zoom, col, row)`` +* ``pixel_x_size``: Horizontal pixel size in SRID units. +* ``pixel_y_size``: Vertical pixel size in SRID units. +* ``left``: Left coordinate. +* ``top``: Top coordinate. +* ``right``: Right coordinate. +* ``bottom``: Bottom coordinate. +* ``width``: Horizontal size in pixels. +* ``height``: Vertical size in pixels. + +### Methods + +#### Get other Tiles +* ``get_parent()``: Returns parent Tile. +* ``get_children()``: Returns children Tiles. +* ``get_neighbors(connectedness=8)``: Returns a maximum of 8 valid neighbor Tiles. + * ``connectedness``: ``4`` or ``8``. Direct neighbors (up to 4) or corner neighbors (up to 8). +* ``intersecting(TilePyramid)``: Return all tiles from other TilePyramid intersecting with tile. This helps translating between TilePyramids with different metatiling +settings. + * ``tilepyramid``: ``TilePyramid`` object. + + +#### Geometry +* ``bounds(pixelbuffer=0)``: Returns a tuple of ``(left, bottom, right, top)`` coordinates. + * ``pixelbuffer``: Optional buffer around tile boundaries in pixels. +* ``bbox(pixelbuffer=0)``: Returns a ``Polygon`` geometry. + * ``pixelbuffer``: Optional buffer around tile boundaries in pixels. +* ``affine(pixelbuffer=0)``: Returns an [``Affine``](https://github.com/sgillies/affine) object. + * ``pixelbuffer``: Optional buffer around tile boundaries in pixels. + + + +#### Other +* ``shape(pixelbuffer=0)``: Returns a tuple with ``(height, width)``. + * ``pixelbuffer``: Optional buffer around tile boundaries in pixels. +* ``is_valid()``: Returns ``True`` if tile ID is valid in this tile matrix and ``False`` if it isn't. diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..1fea590 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,70 @@ +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[project] +name = "tilematrix" +dynamic = ["version"] +description = "helps handling tile pyramids" +readme = "README.rst" +license = "MIT" +authors = [ + { name = "Joachim Ungar", email = "joachim.ungar@gmail.com" }, +] +classifiers = [ + "Development Status :: 5 - Production/Stable", + "Intended Audience :: Developers", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Topic :: Scientific/Engineering :: GIS", +] +dependencies = [ + "affine", + "click", + "geojson", + "rasterio>=1.0.21", + "shapely", +] + +[project.scripts] +tmx = "tilematrix.tmx.main:tmx" + +[project.urls] +Homepage = "https://github.com/ungarj/tilematrix" + +[project.optional-dependencies] +test = [ + "black", + "coveralls", + "flake8", + "pre-commit", + "pytest", + "pytest-cov" +] + +[tool.hatch.version] +path = "tilematrix/__init__.py" + +[tool.hatch.build.targets.sdist] +include = [ + "/tilematrix", +] + +[tool.black] +include = '\.pyi?$' +exclude = ''' +/( + \.git + | \.hg + | \.mypy_cache + | \.tox + | \.venv + | _build + | buck-out + | build + | dist +)/ +''' diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 0000000..3c97111 --- /dev/null +++ b/pytest.ini @@ -0,0 +1,2 @@ +[pytest] +addopts = --durations 20 --verbose --nf --cov=tilematrix \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..6988760 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,6 @@ +affine +click>=7.1.1 +geojson +rasterio>=1.0.21 +shapely +tilematrix \ No newline at end of file diff --git a/test/conftest.py b/test/conftest.py new file mode 100644 index 0000000..eaf8f46 --- /dev/null +++ b/test/conftest.py @@ -0,0 +1,299 @@ +"""Fixtures for test suite.""" + +import pytest +from shapely.geometry import Point, Polygon, shape + +example_proj = """ + +proj=ortho + +lat_0=-90 + +lon_0=0 + +x_0=0 + +y_0=0 + +ellps=WGS84 + +units=m +no_defs + """ + + +@pytest.fixture +def grid_definition_proj(): + """Custom grid definition using a proj string.""" + return { + "shape": (1, 1), + "bounds": (-4000000.0, -4000000.0, 4000000.0, 4000000.0), + "is_global": False, + "proj": example_proj, + } + + +@pytest.fixture +def grid_definition_epsg(): + """Custom grid definition using an EPSG code.""" + return { + "shape": (1, 1), + "bounds": (2426378.0132, 1528101.2618, 6293974.6215, 5395697.8701), + "is_global": False, + "epsg": 3035, + } + + +@pytest.fixture +def proj(): + return example_proj + + +@pytest.fixture +def wkt(): + return """ + PROJCS["ETRS89 / LAEA Europe", + GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989", + SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]], + TOWGS84[0,0,0,0,0,0,0], + AUTHORITY["EPSG","6258"]], + PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], + UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]], + AUTHORITY["EPSG","4258"]], + PROJECTION["Lambert_Azimuthal_Equal_Area"], + PARAMETER["latitude_of_center",52], + PARAMETER["longitude_of_center",10], + PARAMETER["false_easting",4321000], + PARAMETER["false_northing",3210000], + UNIT["metre",1, + AUTHORITY["EPSG","9001"]], + AUTHORITY["EPSG","3035"]] + """ + + +@pytest.fixture +def invalid_geom(): + return Polygon( + [ + (0, 0), + (0, 3), + (3, 3), + (3, 0), + (2, 0), + (2, 2), + (1, 2), + (1, 1), + (2, 1), + (2, 0), + (0, 0), + ] + ) + + +@pytest.fixture +def grid_definition_irregular(): + return { + "shape": [161, 315], + "bounds": [141920, 89840, 948320, 502000], + "is_global": False, + "srs": {"epsg": 31259}, + } + + +@pytest.fixture +def tile_bbox(): + return Polygon( + [ + [-163.125, 73.125], + [-157.5, 73.125], + [-157.5, 67.5], + [-163.125, 67.5], + [-163.125, 73.125], + ] + ) + + +@pytest.fixture +def tile_bbox_buffer(): + return Polygon( + [ + [-163.14697265625, 73.14697265625], + [-157.47802734375, 73.14697265625], + [-157.47802734375, 67.47802734375], + [-163.14697265625, 67.47802734375], + [-163.14697265625, 73.14697265625], + ] + ) + + +@pytest.fixture +def tl_polygon(): + return shape( + { + "type": "Polygon", + "coordinates": ( + ( + (-174.35302734375, 84.35302734375), + (-174.35302734375, 90.0), + (-180.02197265625, 90.0), + (-180.02197265625, 84.35302734375), + (-174.35302734375, 84.35302734375), + ), + ), + } + ) + + +@pytest.fixture +def bl_polygon(): + return shape( + { + "type": "Polygon", + "coordinates": ( + ( + (-174.35302734375, -90.0), + (-174.35302734375, -84.35302734375), + (-180.02197265625, -84.35302734375), + (-180.02197265625, -90.0), + (-174.35302734375, -90.0), + ), + ), + } + ) + + +@pytest.fixture +def overflow_polygon(): + return shape( + { + "type": "Polygon", + "coordinates": ( + ( + (0.703125, -90.0), + (0.703125, 90.0), + (-180.703125, 90.0), + (-180.703125, -90.0), + (0.703125, -90.0), + ), + ), + } + ) + + +@pytest.fixture +def point(): + return Point(16.36, 48.20) + + +@pytest.fixture +def multipoint(): + return shape( + { + "type": "MultiPoint", + "coordinates": [ + (14.464033917048539, 50.08528287347832), + (16.364693096743736, 48.20196113681686), + ], + } + ) + + +@pytest.fixture +def linestring(): + return shape( + { + "type": "LineString", + "coordinates": [ + (8.219788038779399, 48.04680919045518), + (8.553359409223447, 47.98081838641845), + (9.41408206547689, 48.13835399026023), + (10.71989383306024, 48.64871043557477), + (11.683555942439085, 48.794127916044104), + (12.032991977596737, 49.02749868427421), + ], + } + ) + + +@pytest.fixture +def multilinestring(): + return shape( + { + "type": "MultiLineString", + "coordinates": [ + [ + (8.219788038779399, 48.04680919045518), + (8.553359409223447, 47.98081838641845), + (9.41408206547689, 48.13835399026023), + (10.71989383306024, 48.64871043557477), + (11.683555942439085, 48.794127916044104), + (12.032991977596737, 49.02749868427421), + ], + [ + (33.206893344868945, 0.261534735511418), + (33.18725630059802, 0.428191229652711), + (32.8931140479927, 1.31144481038541), + (32.80150465264725, 1.366544806316611), + (32.62475833510098, 1.471712805584616), + (32.51003665541302, 1.536754055177965), + (32.36248752211165, 1.606878973798047), + ], + ], + } + ) + + +@pytest.fixture +def polygon(): + return shape( + { + "type": "Polygon", + "coordinates": [ + [ + (8.219788038779399, 48.04680919045518), + (8.553359409223447, 47.98081838641845), + (9.41408206547689, 48.13835399026023), + (10.71989383306024, 48.64871043557477), + (11.683555942439085, 48.794127916044104), + (12.032991977596737, 49.02749868427421), + (8.219788038779399, 48.04680919045518), + ] + ], + } + ) + + +@pytest.fixture +def multipolygon(): + return shape( + { + "type": "MultiPolygon", + "coordinates": [ + [ + [ + (8.219788038779399, 48.04680919045518), + (8.553359409223447, 47.98081838641845), + (9.41408206547689, 48.13835399026023), + (10.71989383306024, 48.64871043557477), + (11.683555942439085, 48.794127916044104), + (12.032991977596737, 49.02749868427421), + ] + ], + [ + [ + (33.206893344868945, 0.261534735511418), + (33.18725630059802, 0.428191229652711), + (32.8931140479927, 1.31144481038541), + (32.80150465264725, 1.366544806316611), + (32.62475833510098, 1.471712805584616), + (32.51003665541302, 1.536754055177965), + (32.36248752211165, 1.606878973798047), + ] + ], + ], + } + ) + + +@pytest.fixture +def tile_bounds_polygon(): + return shape( + { + "type": "Polygon", + "coordinates": [ + [(0, 0), (0, 45), (22.5, 45), (22.5, 22.5), (45, 22.5), (45, 0), (0, 0)] + ], + } + ) diff --git a/test/requirements.txt b/test/requirements.txt new file mode 100644 index 0000000..3b573da --- /dev/null +++ b/test/requirements.txt @@ -0,0 +1,6 @@ +black +coveralls +flake8 +pre-commit +pytest +pytest-cov diff --git a/test/test_cli.py b/test/test_cli.py new file mode 100644 index 0000000..7b9d6c7 --- /dev/null +++ b/test/test_cli.py @@ -0,0 +1,262 @@ +import geojson +from click.testing import CliRunner +from shapely import wkt +from shapely.geometry import shape + +from tilematrix import TilePyramid, __version__ +from tilematrix.tmx.main import tmx + + +def test_version(): + result = CliRunner().invoke(tmx, ["--version"]) + assert result.exit_code == 0 + assert __version__ in result.output + + +def test_bounds(): + zoom, row, col = 12, 0, 0 + for grid in ["geodetic", "mercator"]: + _run_bbox_bounds(zoom, row, col, "bounds", grid=grid) + for metatiling in [1, 2, 4, 8, 16]: + _run_bbox_bounds(zoom, row, col, "bounds", metatiling=metatiling) + for pixelbuffer in [0, 1, 10]: + _run_bbox_bounds(zoom, row, col, "bounds", pixelbuffer=pixelbuffer) + for tile_size in [256, 512, 1024]: + _run_bbox_bounds(zoom, row, col, "bounds", tile_size=tile_size) + + +def test_bbox(): + zoom, row, col = 12, 0, 0 + for grid in ["geodetic", "mercator"]: + _run_bbox_bounds(zoom, row, col, "bbox", grid=grid) + for metatiling in [1, 2, 4, 8, 16]: + _run_bbox_bounds(zoom, row, col, "bbox", metatiling=metatiling) + for pixelbuffer in [0, 1, 10]: + _run_bbox_bounds(zoom, row, col, "bbox", pixelbuffer=pixelbuffer) + for tile_size in [256, 512, 1024]: + _run_bbox_bounds(zoom, row, col, "bbox", tile_size=tile_size) + for output_format in ["WKT", "GeoJSON"]: + _run_bbox_bounds(zoom, row, col, "bbox", output_format=output_format) + + +def test_tile(): + point = [10, 20] + zoom = 5 + for grid in ["geodetic", "mercator"]: + _run_tile(zoom, point, grid=grid) + for metatiling in [1, 2, 4, 8, 16]: + _run_tile(zoom, point, metatiling=metatiling) + for pixelbuffer in [0, 1, 10]: + _run_tile(zoom, point, pixelbuffer=pixelbuffer) + for tile_size in [256, 512, 1024]: + _run_tile(zoom, point, tile_size=tile_size) + for output_format in ["Tile", "WKT", "GeoJSON"]: + _run_tile(zoom, point, output_format=output_format) + + +def test_tiles(): + bounds = (10, 20, 30, 40) + zoom = 6 + for grid in ["geodetic", "mercator"]: + _run_tiles(zoom, bounds, grid=grid) + for metatiling in [1, 2, 4, 8, 16]: + _run_tiles(zoom, bounds, metatiling=metatiling) + for pixelbuffer in [0, 1, 10]: + _run_tiles(zoom, bounds, pixelbuffer=pixelbuffer) + for tile_size in [256, 512, 1024]: + _run_tiles(zoom, bounds, tile_size=tile_size) + for output_format in ["Tile", "WKT", "GeoJSON"]: + _run_tiles(zoom, bounds, output_format=output_format) + + +def test_snap_bounds(): + zoom = 6 + bounds = (10, 20, 30, 40) + for grid in ["geodetic", "mercator"]: + _run_snap_bounds(zoom, bounds, "snap-bounds", grid=grid) + for metatiling in [1, 2, 4, 8, 16]: + _run_snap_bounds(zoom, bounds, "snap-bounds", metatiling=metatiling) + for pixelbuffer in [0, 1, 10]: + _run_snap_bounds(zoom, bounds, "snap-bounds", pixelbuffer=pixelbuffer) + for tile_size in [256, 512, 1024]: + _run_snap_bounds(zoom, bounds, "snap-bounds", tile_size=tile_size) + + +def test_snap_bbox(): + zoom = 6 + bounds = (10, 20, 30, 40) + for grid in ["geodetic", "mercator"]: + _run_snap_bounds(zoom, bounds, "snap-bbox", grid=grid) + for metatiling in [1, 2, 4, 8, 16]: + _run_snap_bounds(zoom, bounds, "snap-bbox", metatiling=metatiling) + for pixelbuffer in [0, 1, 10]: + _run_snap_bounds(zoom, bounds, "snap-bbox", pixelbuffer=pixelbuffer) + for tile_size in [256, 512, 1024]: + _run_snap_bounds(zoom, bounds, "snap-bbox", tile_size=tile_size) + + +def _run_bbox_bounds( + zoom, + row, + col, + command=None, + grid="geodetic", + metatiling=1, + pixelbuffer=0, + tile_size=256, + output_format="WKT", +): + tile = TilePyramid(grid, metatiling=metatiling, tile_size=tile_size).tile( + zoom, row, col + ) + result = CliRunner().invoke( + tmx, + [ + "--pixelbuffer", + str(pixelbuffer), + "--metatiling", + str(metatiling), + "--grid", + grid, + "--tile_size", + str(tile_size), + "--output_format", + output_format, + command, + str(zoom), + str(row), + str(col), + ], + ) + assert result.exit_code == 0 + if command == "bounds": + assert result.output.strip() == " ".join(map(str, tile.bounds(pixelbuffer))) + elif output_format == "WKT": + assert wkt.loads(result.output.strip()).almost_equals(tile.bbox(pixelbuffer)) + elif output_format == "GeoJSON": + assert shape(geojson.loads(result.output.strip())).almost_equals( + tile.bbox(pixelbuffer) + ) + + +def _run_tile( + zoom, + point, + grid="geodetic", + metatiling=1, + pixelbuffer=0, + tile_size=256, + output_format="WKT", +): + ( + x, + y, + ) = point + tile = TilePyramid(grid, metatiling=metatiling, tile_size=tile_size).tile_from_xy( + x, y, zoom + ) + result = CliRunner().invoke( + tmx, + [ + "--pixelbuffer", + str(pixelbuffer), + "--metatiling", + str(metatiling), + "--grid", + grid, + "--tile_size", + str(tile_size), + "--output_format", + output_format, + "tile", + str(zoom), + str(x), + str(y), + ], + ) + assert result.exit_code == 0 + if output_format == "Tile": + assert result.output.strip() == " ".join(map(str, tile.id)) + elif output_format == "WKT": + assert wkt.loads(result.output.strip()).almost_equals(tile.bbox(pixelbuffer)) + elif output_format == "GeoJSON": + feature = geojson.loads(result.output.strip())["features"][0] + assert shape(feature["geometry"]).almost_equals(tile.bbox(pixelbuffer)) + + +def _run_tiles( + zoom, + bounds, + grid="geodetic", + metatiling=1, + pixelbuffer=0, + tile_size=256, + output_format="WKT", +): + left, bottom, right, top = bounds + tiles = list( + TilePyramid(grid, metatiling=metatiling, tile_size=tile_size).tiles_from_bounds( + bounds, zoom + ) + ) + result = CliRunner().invoke( + tmx, + [ + "--pixelbuffer", + str(pixelbuffer), + "--metatiling", + str(metatiling), + "--grid", + grid, + "--tile_size", + str(tile_size), + "--output_format", + output_format, + "tiles", + str(zoom), + str(left), + str(bottom), + str(right), + str(top), + ], + ) + assert result.exit_code == 0 + if output_format == "Tile": + assert result.output.count("\n") == len(tiles) + elif output_format == "WKT": + assert result.output.count("\n") == len(tiles) + elif output_format == "GeoJSON": + features = geojson.loads(result.output.strip())["features"] + assert len(features) == len(tiles) + + +def _run_snap_bounds( + zoom, + bounds, + command=None, + grid="geodetic", + metatiling=1, + pixelbuffer=0, + tile_size=256, +): + left, bottom, right, top = bounds + result = CliRunner().invoke( + tmx, + [ + "--pixelbuffer", + str(pixelbuffer), + "--metatiling", + str(metatiling), + "--grid", + grid, + "--tile_size", + str(tile_size), + command, + str(zoom), + str(left), + str(bottom), + str(right), + str(top), + ], + ) + assert result.exit_code == 0 diff --git a/test/test_dump_load.py b/test/test_dump_load.py new file mode 100644 index 0000000..feaad54 --- /dev/null +++ b/test/test_dump_load.py @@ -0,0 +1,26 @@ +from tilematrix import TilePyramid + + +def test_geodetic(): + tp = TilePyramid("geodetic", metatiling=2) + tp_dict = tp.to_dict() + assert isinstance(tp_dict, dict) + tp2 = TilePyramid.from_dict(tp_dict) + assert tp == tp2 + + +def test_mercator(): + tp = TilePyramid("mercator", metatiling=4) + tp_dict = tp.to_dict() + assert isinstance(tp_dict, dict) + tp2 = TilePyramid.from_dict(tp_dict) + assert tp == tp2 + + +def test_custom(grid_definition_proj, grid_definition_epsg): + for grid_def in [grid_definition_proj, grid_definition_epsg]: + tp = TilePyramid(grid_def, metatiling=8) + tp_dict = tp.to_dict() + assert isinstance(tp_dict, dict) + tp2 = TilePyramid.from_dict(tp_dict) + assert tp == tp2 diff --git a/test/test_geometries.py b/test/test_geometries.py new file mode 100644 index 0000000..3c8fb9d --- /dev/null +++ b/test/test_geometries.py @@ -0,0 +1,432 @@ +"""Tile geometries and tiles from geometries.""" + +from types import GeneratorType + +import pytest +from shapely.geometry import Point, Polygon, shape + +from tilematrix import Tile, TilePyramid + + +def test_top_left_coord(): + """Top left coordinate.""" + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert (tile.left, tile.top) == (-163.125, 73.125) + + +def test_tile_bbox(tile_bbox): + """Tile bounding box.""" + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bbox().equals(tile_bbox) + + +def test_tile_bbox_buffer(tile_bbox_buffer, tl_polygon, bl_polygon, overflow_polygon): + """Tile bounding box with buffer.""" + # default + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bbox(1).equals(tile_bbox_buffer) + + # first row of tilematrix + tile = tp.tile(5, 0, 0) + assert tile.bbox(1) == tl_polygon + + # last row of tilematrix + tile = tp.tile(5, 31, 0) + assert tile.bbox(1) == bl_polygon + + # overflowing all tilepyramid bounds + tile = tp.tile(0, 0, 0) + assert tile.bbox(1) == overflow_polygon + + +def test_tile_bounds(): + """Tile bounds.""" + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bounds() == (-163.125, 67.5, -157.5, 73.125) + + +def test_tile_bounds_buffer(): + """Tile bounds with buffer.""" + tp = TilePyramid("geodetic") + # default + tile = tp.tile(5, 3, 3) + testbounds = (-163.14697265625, 67.47802734375, -157.47802734375, 73.14697265625) + assert tile.bounds(1) == testbounds + + # first row of tilematrix + tile = tp.tile(5, 0, 0) + testbounds = (-180.02197265625, 84.35302734375, -174.35302734375, 90.0) + assert tile.bounds(1) == testbounds + + # last row of tilematrix + tile = tp.tile(5, 31, 0) + testbounds = (-180.02197265625, -90.0, -174.35302734375, -84.35302734375) + assert tile.bounds(1) == testbounds + + # overflowing all tilepyramid bounds + tile = tp.tile(0, 0, 0) + testbounds = (-180.703125, -90.0, 0.703125, 90.0) + assert tile.bounds(1) == testbounds + + +def test_tiles_from_bounds(): + """Get tiles intersecting with given bounds.""" + tp = TilePyramid("geodetic") + # valid bounds + bounds = (-163.125, 67.5, -157.5, 73.125) + control_tiles = {(5, 3, 3)} + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 5)} + assert control_tiles == test_tiles + # invalid bounds + try: + {tile.id for tile in tp.tiles_from_bounds((3, 5), 5)} + raise Exception() + except ValueError: + pass + # cross the antimeridian on the western side + bounds = (-183.125, 67.5, -177.5, 73.125) + control_tiles = {(5, 3, 0), (5, 3, 63)} + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 5)} + assert control_tiles == test_tiles + # cross the antimeridian on the eastern side + bounds = (177.5, 67.5, 183.125, 73.125) + control_tiles = {(5, 3, 0), (5, 3, 63)} + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 5)} + assert control_tiles == test_tiles + # cross the antimeridian on both sudes + bounds = (-183, 67.5, 183.125, 73.125) + control_tiles = { + (3, 0, 0), + (3, 0, 1), + (3, 0, 2), + (3, 0, 3), + (3, 0, 4), + (3, 0, 5), + (3, 0, 6), + (3, 0, 7), + (3, 0, 8), + (3, 0, 9), + (3, 0, 10), + (3, 0, 11), + (3, 0, 12), + (3, 0, 13), + (3, 0, 14), + (3, 0, 15), + } + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 3)} + assert control_tiles == test_tiles + + +def test_tiles_from_bbox(): + """Get tiles intersecting with bounding box.""" + test_bbox = shape( + { + "type": "Polygon", + "coordinates": [ + [ + (5.625, 61.875), + (56.25, 61.875), + (56.25, 28.125), + (5.625, 28.125), + (5.625, 28.125), + (5.625, 61.875), + ] + ], + } + ) + test_tiles = { + (5, 5, 33), + (5, 6, 33), + (5, 7, 33), + (5, 8, 33), + (5, 9, 33), + (5, 10, 33), + (5, 5, 34), + (5, 6, 34), + (5, 7, 34), + (5, 8, 34), + (5, 9, 34), + (5, 10, 34), + (5, 5, 35), + (5, 6, 35), + (5, 7, 35), + (5, 8, 35), + (5, 9, 35), + (5, 10, 35), + (5, 5, 36), + (5, 6, 36), + (5, 7, 36), + (5, 8, 36), + (5, 9, 36), + (5, 10, 36), + (5, 5, 37), + (5, 6, 37), + (5, 7, 37), + (5, 8, 37), + (5, 9, 37), + (5, 10, 37), + (5, 5, 38), + (5, 6, 38), + (5, 7, 38), + (5, 8, 38), + (5, 9, 38), + (5, 10, 38), + (5, 5, 39), + (5, 6, 39), + (5, 7, 39), + (5, 8, 39), + (5, 9, 39), + (5, 10, 39), + (5, 5, 40), + (5, 6, 40), + (5, 7, 40), + (5, 8, 40), + (5, 9, 40), + (5, 10, 40), + (5, 5, 41), + (5, 6, 41), + (5, 7, 41), + (5, 8, 41), + (5, 9, 41), + (5, 10, 41), + } + tp = TilePyramid("geodetic") + bbox_tiles = {tile.id for tile in tp.tiles_from_bbox(test_bbox, 5)} + assert test_tiles == bbox_tiles + + +def test_tiles_from_empty_geom(): + """Get tiles from empty geometry.""" + test_geom = Polygon() + tp = TilePyramid("geodetic") + empty_tiles = {tile.id for tile in tp.tiles_from_geom(test_geom, 6)} + assert empty_tiles == set([]) + + +def test_tiles_from_invalid_geom(invalid_geom): + """Get tiles from empty geometry.""" + tp = TilePyramid("geodetic") + with pytest.raises(ValueError): + list(tp.tiles_from_geom(invalid_geom, 6)) + + +def test_tiles_from_point(point): + """Get tile from point.""" + for metatiling in [1, 2, 4, 8, 16]: + tp = TilePyramid("geodetic", metatiling=metatiling) + tile_bbox = next(tp.tiles_from_geom(point, 6)).bbox() + assert point.within(tile_bbox) + tp = TilePyramid("geodetic") + with pytest.raises(ValueError): + next(tp.tiles_from_geom(Point(-300, 100), 6)) + + +def test_tiles_from_multipoint(multipoint): + """Get tiles from multiple points.""" + test_tiles = {(9, 113, 553), (9, 118, 558)} + tp = TilePyramid("geodetic") + multipoint_tiles = {tile.id for tile in tp.tiles_from_geom(multipoint, 9)} + assert multipoint_tiles == test_tiles + + +def test_tiles_from_linestring(linestring): + """Get tiles from LineString.""" + test_tiles = { + (8, 58, 270), + (8, 58, 271), + (8, 58, 272), + (8, 58, 273), + (8, 59, 267), + (8, 59, 268), + (8, 59, 269), + (8, 59, 270), + } + tp = TilePyramid("geodetic") + linestring_tiles = {tile.id for tile in tp.tiles_from_geom(linestring, 8)} + assert linestring_tiles == test_tiles + + +def test_tiles_from_multilinestring(multilinestring): + """Get tiles from MultiLineString.""" + test_tiles = { + (8, 58, 270), + (8, 58, 271), + (8, 58, 272), + (8, 58, 273), + (8, 59, 267), + (8, 59, 268), + (8, 59, 269), + (8, 59, 270), + (8, 125, 302), + (8, 126, 302), + (8, 126, 303), + (8, 127, 303), + } + tp = TilePyramid("geodetic") + multilinestring_tiles = {tile.id for tile in tp.tiles_from_geom(multilinestring, 8)} + assert multilinestring_tiles == test_tiles + + +def test_tiles_from_polygon(polygon): + """Get tiles from Polygon.""" + test_tiles = { + (9, 116, 544), + (9, 116, 545), + (9, 116, 546), + (9, 117, 540), + (9, 117, 541), + (9, 117, 542), + (9, 117, 543), + (9, 117, 544), + (9, 117, 545), + (9, 118, 536), + (9, 118, 537), + (9, 118, 538), + (9, 118, 539), + (9, 118, 540), + (9, 118, 541), + (9, 119, 535), + (9, 119, 536), + (9, 119, 537), + (9, 119, 538), + } + tp = TilePyramid("geodetic") + polygon_tiles = {tile.id for tile in tp.tiles_from_geom(polygon, 9)} + assert polygon_tiles == test_tiles + + +def test_tiles_from_multipolygon(multipolygon): + """Get tiles from MultiPolygon.""" + test_tiles = { + (9, 116, 544), + (9, 116, 545), + (9, 116, 546), + (9, 117, 540), + (9, 117, 541), + (9, 117, 542), + (9, 117, 543), + (9, 117, 544), + (9, 117, 545), + (9, 118, 536), + (9, 118, 537), + (9, 118, 538), + (9, 118, 539), + (9, 118, 540), + (9, 118, 541), + (9, 119, 535), + (9, 119, 536), + (9, 119, 537), + (9, 119, 538), + (9, 251, 604), + (9, 251, 605), + (9, 252, 604), + (9, 252, 605), + (9, 253, 605), + (9, 253, 606), + (9, 254, 605), + (9, 254, 606), + (9, 255, 606), + } + tp = TilePyramid("geodetic") + multipolygon_tiles = {tile.id for tile in tp.tiles_from_geom(multipolygon, 9)} + assert multipolygon_tiles == test_tiles + + +def test_tiles_from_point_batches(point): + """Get tile from point.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(point, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(point, zoom))) + + +def test_tiles_from_multipoint_batches(multipoint): + """Get tiles from multiple points.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(multipoint, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(multipoint, zoom))) + + +def test_tiles_from_linestring_batches(linestring): + """Get tiles from LineString.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(linestring, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(linestring, zoom))) + + +def test_tiles_from_multilinestring_batches(multilinestring): + """Get tiles from MultiLineString.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(multilinestring, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(multilinestring, zoom))) + + +def test_tiles_from_polygon_batches(polygon): + """Get tiles from Polygon.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(polygon, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(polygon, zoom))) + + +def test_tiles_from_multipolygon_batches(multipolygon): + """Get tiles from MultiPolygon.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(multipolygon, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(multipolygon, zoom))) diff --git a/test/test_grids.py b/test/test_grids.py new file mode 100644 index 0000000..6d09f9c --- /dev/null +++ b/test/test_grids.py @@ -0,0 +1,108 @@ +"""TilePyramid creation.""" + +import math + +import pytest +from shapely.geometry import box + +from tilematrix import PYRAMID_PARAMS, GridDefinition, TilePyramid + + +def test_grid_init(grid_definition_proj): + grid = GridDefinition(**dict(PYRAMID_PARAMS["geodetic"], grid="custom")) + custom_grid = TilePyramid(grid_definition_proj).grid + # make sure standard grid gets detected + assert grid.type == "geodetic" + # create grid from grid + assert GridDefinition(grid) + # create grid from dict + assert GridDefinition.from_dict(grid.to_dict()) + # __repr__ + assert str(grid) + assert str(custom_grid) + # __hash__ + assert hash(grid) + assert hash(custom_grid) + + +def test_deprecated(): + grid = TilePyramid("geodetic").grid + assert grid.srid + + +def test_proj_str(grid_definition_proj): + """Initialize with proj string.""" + tp = TilePyramid(grid_definition_proj) + assert tp.tile(0, 0, 0).bounds() == grid_definition_proj["bounds"] + + +def test_epsg_code(grid_definition_epsg): + """Initialize with EPSG code.""" + tp = TilePyramid(grid_definition_epsg) + assert tp.tile(0, 0, 0).bounds() == grid_definition_epsg["bounds"] + + +def test_shape_error(grid_definition_epsg): + """Raise error when shape aspect ratio is not bounds apsect ratio.""" + grid_definition_epsg.update( + bounds=(2426378.0132, 1528101.2618, 6293974.6215, 5446513.5222) + ) + with pytest.raises(ValueError): + TilePyramid(grid_definition_epsg) + + +def test_neighbors(grid_definition_epsg): + """Initialize with EPSG code.""" + tp = TilePyramid(grid_definition_epsg) + neighbor_ids = set([t.id for t in tp.tile(1, 0, 0).get_neighbors()]) + control_ids = set([(1, 1, 0), (1, 0, 1), (1, 1, 1)]) + assert len(neighbor_ids.symmetric_difference(control_ids)) == 0 + + +def test_irregular_grids(grid_definition_irregular): + for metatiling in [1, 2, 4, 8]: + for pixelbuffer in [0, 100]: + tp = TilePyramid(grid_definition_irregular, metatiling=metatiling) + assert tp.matrix_height(0) == math.ceil(161 / metatiling) + assert tp.matrix_width(0) == math.ceil(315 / metatiling) + assert tp.pixel_x_size(0) == tp.pixel_y_size(0) + for tile in [ + tp.tile(0, 0, 0), + tp.tile(0, 10, 0), + tp.tile(0, tp.matrix_height(0) - 1, tp.matrix_width(0) - 1), + ]: + # pixel sizes have to be squares + assert tile.pixel_x_size == tile.pixel_y_size + assert tile.pixel_x_size == 10.0 + # pixelbuffer yields different shape and bounds + assert tile.shape(10) != tile.shape() + assert tile.bounds(10) != tile.bounds() + # without pixelbuffers, tile bounds have to be inside TilePyramid bounds + assert tile.left >= tp.left + assert tile.bottom >= tp.bottom + assert tile.right <= tp.right + assert tile.top <= tp.top + + # with pixelbuffers, some tile bounds have to be outside TilePyramid bounds + if pixelbuffer: + # tile on top left corner + tile = tp.tile(0, 0, 0) + assert tile.bounds(pixelbuffer).left < tp.left + assert tile.bounds(pixelbuffer).top > tp.top + + # tile on lower right corner + tile = tp.tile(0, tp.matrix_height(0) - 1, tp.matrix_width(0) - 1) + print(tp) + assert tile.bounds(pixelbuffer).bottom < tp.bottom + assert tile.bounds(pixelbuffer).right > tp.right + + +def test_tiles_from_bounds(grid_definition_irregular): + bounds = (755336.179, 300068.615, 791558.022, 319499.955) + bbox = box(*bounds) + tp = TilePyramid(grid_definition_irregular, metatiling=4) + tiles_bounds = list(tp.tiles_from_bounds(bounds, 0)) + tiles_bbox = list(tp.tiles_from_bbox(bbox, 0)) + tiles_geom = list(tp.tiles_from_geom(bbox, 0)) + + assert set(tiles_bounds) == set(tiles_bbox) == set(tiles_geom) diff --git a/test/test_helper_funcs.py b/test/test_helper_funcs.py new file mode 100644 index 0000000..92d7fc5 --- /dev/null +++ b/test/test_helper_funcs.py @@ -0,0 +1,99 @@ +import pytest +from rasterio.crs import CRS +from shapely.geometry import box + +from tilematrix import TilePyramid, clip_geometry_to_srs_bounds, validate_zoom +from tilematrix._funcs import _get_crs, _verify_shape_bounds + + +def test_antimeridian_clip(invalid_geom): + """Clip on antimeridian.""" + tp = TilePyramid("geodetic") + tp_bounds = box(tp.left, tp.bottom, tp.right, tp.top) + + # extends on the western side + geometry = box(-183.125, 67.5, -177.5, 73.125) + # get GeometryCollection + out_geom = clip_geometry_to_srs_bounds(geometry, tp) + assert out_geom.geom_type == "GeometryCollection" + # get list + out_geom = clip_geometry_to_srs_bounds(geometry, tp, multipart=True) + assert isinstance(out_geom, list) + for sub_geom in out_geom: + assert sub_geom.within(tp_bounds) + + # extends on the eastern side + geometry = box(177.5, 67.5, 183.125, 73.125) + # get GeometryCollection + out_geom = clip_geometry_to_srs_bounds(geometry, tp) + assert out_geom.geom_type == "GeometryCollection" + # get list + out_geom = clip_geometry_to_srs_bounds(geometry, tp, multipart=True) + assert isinstance(out_geom, list) + for sub_geom in out_geom: + assert sub_geom.within(tp_bounds) + + # extends on both sides + geometry = box(-183.125, 67.5, 183.125, 73.125) + # get GeometryCollection + out_geom = clip_geometry_to_srs_bounds(geometry, tp) + assert out_geom.geom_type == "GeometryCollection" + # get list + out_geom = clip_geometry_to_srs_bounds(geometry, tp, multipart=True) + assert isinstance(out_geom, list) + for sub_geom in out_geom: + assert sub_geom.within(tp_bounds) + assert len(out_geom) == 3 + + # fail on invalid geometry + with pytest.raises(ValueError): + clip_geometry_to_srs_bounds(invalid_geom, tp) + + +def test_no_clip(): + """Geometry is within TilePyramid bounds.""" + tp = TilePyramid("geodetic") + geometry = box(177.5, 67.5, -177.5, 73.125) + + # no multipart + out_geom = clip_geometry_to_srs_bounds(geometry, tp) + assert geometry == out_geom + + # multipart + out_geom = clip_geometry_to_srs_bounds(geometry, tp, multipart=True) + assert isinstance(out_geom, list) + assert len(out_geom) == 1 + assert geometry == out_geom[0] + + +def test_validate_zoom(): + with pytest.raises(TypeError): + validate_zoom(5.0) + with pytest.raises(ValueError): + validate_zoom(-1) + + +def test_verify_shape_bounds(): + # invalid shape + with pytest.raises(TypeError): + _verify_shape_bounds((1, 2, 3), (1, 2, 3, 4)) + # invalid bounds + with pytest.raises(TypeError): + _verify_shape_bounds((1, 2), (1, 2, 3, 4, 5)) + + _verify_shape_bounds((2, 1), (1, 2, 3, 6)) + + +def test_get_crs(proj, wkt): + # no dictionary + with pytest.raises(TypeError): + _get_crs("no dict") + # valid WKT + assert isinstance(_get_crs(dict(wkt=wkt)), CRS) + # valid EPSG + assert isinstance(_get_crs(dict(epsg=3857)), CRS) + # valid PROJ + assert isinstance(_get_crs(dict(proj=proj)), CRS) + # none of above + with pytest.raises(TypeError): + _get_crs(dict(something_else=None)) diff --git a/test/test_matrix_shapes.py b/test/test_matrix_shapes.py new file mode 100644 index 0000000..93ff836 --- /dev/null +++ b/test/test_matrix_shapes.py @@ -0,0 +1,150 @@ +#!/usr/bin/env python +"""Test tile matrix shapes.""" + +from tilematrix import TilePyramid + + +def test_geodetic_matrix_shapes(): + """Test shapes of geodetic tile matrices.""" + # 1 metatiling + matrix_shapes = { + 0: (2, 1), + 1: (4, 2), + 2: (8, 4), + 3: (16, 8), + 4: (32, 16), + 5: (64, 32), + 6: (128, 64), + } + tp = TilePyramid("geodetic") + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 2 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (2, 1), + 2: (4, 2), + 3: (8, 4), + 4: (16, 8), + 5: (32, 16), + 6: (64, 32), + } + tp = TilePyramid("geodetic", metatiling=2) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 4 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (2, 1), + 3: (4, 2), + 4: (8, 4), + 5: (16, 8), + 6: (32, 16), + } + tp = TilePyramid("geodetic", metatiling=4) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 8 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (1, 1), + 3: (2, 1), + 4: (4, 2), + 5: (8, 4), + 6: (16, 8), + } + tp = TilePyramid("geodetic", metatiling=8) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 16 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (1, 1), + 3: (1, 1), + 4: (2, 1), + 5: (4, 2), + 6: (8, 4), + } + tp = TilePyramid("geodetic", metatiling=16) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + +def test_mercator_matrix_shapes(): + """Test shapes of mercator tile matrices.""" + # 1 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (2, 2), + 2: (4, 4), + 3: (8, 8), + 4: (16, 16), + 5: (32, 32), + 6: (64, 64), + } + tp = TilePyramid("mercator") + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 2 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (2, 2), + 3: (4, 4), + 4: (8, 8), + 5: (16, 16), + 6: (32, 32), + } + tp = TilePyramid("mercator", metatiling=2) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 4 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (1, 1), + 3: (2, 2), + 4: (4, 4), + 5: (8, 8), + 6: (16, 16), + } + tp = TilePyramid("mercator", metatiling=4) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 8 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (1, 1), + 3: (1, 1), + 4: (2, 2), + 5: (4, 4), + 6: (8, 8), + } + tp = TilePyramid("mercator", metatiling=8) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape + + # 16 metatiling + matrix_shapes = { + 0: (1, 1), + 1: (1, 1), + 2: (1, 1), + 3: (1, 1), + 4: (1, 1), + 5: (2, 2), + 6: (4, 4), + } + tp = TilePyramid("mercator", metatiling=16) + for zoom, shape in matrix_shapes.items(): + assert (tp.matrix_width(zoom), tp.matrix_height(zoom)) == shape diff --git a/test/test_pyramid_cut.py b/test/test_pyramid_cut.py new file mode 100644 index 0000000..54eb5a5 --- /dev/null +++ b/test/test_pyramid_cut.py @@ -0,0 +1,433 @@ +"""Tile geometries and tiles from geometries.""" + +from types import GeneratorType + +import pytest +from shapely.geometry import Point, Polygon, shape + +from tilematrix import Tile, TilePyramid + + +def test_top_left_coord(): + """Top left coordinate.""" + tp = TilePyramid("geodetic") + + tile = tp.tile(5, 3, 3) + assert (tile.left, tile.top) == (-163.125, 73.125) + + +def test_tile_bbox(tile_bbox): + """Tile bounding box.""" + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bbox().equals(tile_bbox) + + +def test_tile_bbox_buffer(tile_bbox_buffer, tl_polygon, bl_polygon, overflow_polygon): + """Tile bounding box with buffer.""" + # default + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bbox(1).equals(tile_bbox_buffer) + + # first row of tilematrix + tile = tp.tile(5, 0, 0) + assert tile.bbox(1) == tl_polygon + + # last row of tilematrix + tile = tp.tile(5, 31, 0) + assert tile.bbox(1) == bl_polygon + + # overflowing all tilepyramid bounds + tile = tp.tile(0, 0, 0) + assert tile.bbox(1) == overflow_polygon + + +def test_tile_bounds(): + """Tile bounds.""" + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bounds() == (-163.125, 67.5, -157.5, 73.125) + + +def test_tile_bounds_buffer(): + """Tile bounds with buffer.""" + tp = TilePyramid("geodetic") + # default + tile = tp.tile(5, 3, 3) + testbounds = (-163.14697265625, 67.47802734375, -157.47802734375, 73.14697265625) + assert tile.bounds(1) == testbounds + + # first row of tilematrix + tile = tp.tile(5, 0, 0) + testbounds = (-180.02197265625, 84.35302734375, -174.35302734375, 90.0) + assert tile.bounds(1) == testbounds + + # last row of tilematrix + tile = tp.tile(5, 31, 0) + testbounds = (-180.02197265625, -90.0, -174.35302734375, -84.35302734375) + assert tile.bounds(1) == testbounds + + # overflowing all tilepyramid bounds + tile = tp.tile(0, 0, 0) + testbounds = (-180.703125, -90.0, 0.703125, 90.0) + assert tile.bounds(1) == testbounds + + +def test_tiles_from_bounds(): + """Get tiles intersecting with given bounds.""" + tp = TilePyramid("geodetic") + # valid bounds + bounds = (-163.125, 67.5, -157.5, 73.125) + control_tiles = {(5, 3, 3)} + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 5)} + assert control_tiles == test_tiles + # invalid bounds + try: + {tile.id for tile in tp.tiles_from_bounds((3, 5), 5)} + raise Exception() + except ValueError: + pass + # cross the antimeridian on the western side + bounds = (-183.125, 67.5, -177.5, 73.125) + control_tiles = {(5, 3, 0), (5, 3, 63)} + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 5)} + assert control_tiles == test_tiles + # cross the antimeridian on the eastern side + bounds = (177.5, 67.5, 183.125, 73.125) + control_tiles = {(5, 3, 0), (5, 3, 63)} + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 5)} + assert control_tiles == test_tiles + # cross the antimeridian on both sudes + bounds = (-183, 67.5, 183.125, 73.125) + control_tiles = { + (3, 0, 0), + (3, 0, 1), + (3, 0, 2), + (3, 0, 3), + (3, 0, 4), + (3, 0, 5), + (3, 0, 6), + (3, 0, 7), + (3, 0, 8), + (3, 0, 9), + (3, 0, 10), + (3, 0, 11), + (3, 0, 12), + (3, 0, 13), + (3, 0, 14), + (3, 0, 15), + } + test_tiles = {tile.id for tile in tp.tiles_from_bounds(bounds, 3)} + assert control_tiles == test_tiles + + +def test_tiles_from_bbox(): + """Get tiles intersecting with bounding box.""" + test_bbox = shape( + { + "type": "Polygon", + "coordinates": [ + [ + (5.625, 61.875), + (56.25, 61.875), + (56.25, 28.125), + (5.625, 28.125), + (5.625, 28.125), + (5.625, 61.875), + ] + ], + } + ) + test_tiles = { + (5, 5, 33), + (5, 6, 33), + (5, 7, 33), + (5, 8, 33), + (5, 9, 33), + (5, 10, 33), + (5, 5, 34), + (5, 6, 34), + (5, 7, 34), + (5, 8, 34), + (5, 9, 34), + (5, 10, 34), + (5, 5, 35), + (5, 6, 35), + (5, 7, 35), + (5, 8, 35), + (5, 9, 35), + (5, 10, 35), + (5, 5, 36), + (5, 6, 36), + (5, 7, 36), + (5, 8, 36), + (5, 9, 36), + (5, 10, 36), + (5, 5, 37), + (5, 6, 37), + (5, 7, 37), + (5, 8, 37), + (5, 9, 37), + (5, 10, 37), + (5, 5, 38), + (5, 6, 38), + (5, 7, 38), + (5, 8, 38), + (5, 9, 38), + (5, 10, 38), + (5, 5, 39), + (5, 6, 39), + (5, 7, 39), + (5, 8, 39), + (5, 9, 39), + (5, 10, 39), + (5, 5, 40), + (5, 6, 40), + (5, 7, 40), + (5, 8, 40), + (5, 9, 40), + (5, 10, 40), + (5, 5, 41), + (5, 6, 41), + (5, 7, 41), + (5, 8, 41), + (5, 9, 41), + (5, 10, 41), + } + tp = TilePyramid("geodetic") + bbox_tiles = {tile.id for tile in tp.tiles_from_bbox(test_bbox, 5)} + assert test_tiles == bbox_tiles + + +def test_tiles_from_empty_geom(): + """Get tiles from empty geometry.""" + test_geom = Polygon() + tp = TilePyramid("geodetic") + empty_tiles = {tile.id for tile in tp.tiles_from_geom(test_geom, 6)} + assert empty_tiles == set([]) + + +def test_tiles_from_invalid_geom(invalid_geom): + """Get tiles from empty geometry.""" + tp = TilePyramid("geodetic") + with pytest.raises(ValueError): + list(tp.tiles_from_geom(invalid_geom, 6)) + + +def test_tiles_from_point(point): + """Get tile from point.""" + for metatiling in [1, 2, 4, 8, 16]: + tp = TilePyramid("geodetic", metatiling=metatiling) + tile_bbox = next(tp.tiles_from_geom(point, 6)).bbox() + assert point.within(tile_bbox) + tp = TilePyramid("geodetic") + with pytest.raises(ValueError): + next(tp.tiles_from_geom(Point(-300, 100), 6)) + + +def test_tiles_from_multipoint(multipoint): + """Get tiles from multiple points.""" + test_tiles = {(9, 113, 553), (9, 118, 558)} + tp = TilePyramid("geodetic") + multipoint_tiles = {tile.id for tile in tp.tiles_from_geom(multipoint, 9)} + assert multipoint_tiles == test_tiles + + +def test_tiles_from_linestring(linestring): + """Get tiles from LineString.""" + test_tiles = { + (8, 58, 270), + (8, 58, 271), + (8, 58, 272), + (8, 58, 273), + (8, 59, 267), + (8, 59, 268), + (8, 59, 269), + (8, 59, 270), + } + tp = TilePyramid("geodetic") + linestring_tiles = {tile.id for tile in tp.tiles_from_geom(linestring, 8)} + assert linestring_tiles == test_tiles + + +def test_tiles_from_multilinestring(multilinestring): + """Get tiles from MultiLineString.""" + test_tiles = { + (8, 58, 270), + (8, 58, 271), + (8, 58, 272), + (8, 58, 273), + (8, 59, 267), + (8, 59, 268), + (8, 59, 269), + (8, 59, 270), + (8, 125, 302), + (8, 126, 302), + (8, 126, 303), + (8, 127, 303), + } + tp = TilePyramid("geodetic") + multilinestring_tiles = {tile.id for tile in tp.tiles_from_geom(multilinestring, 8)} + assert multilinestring_tiles == test_tiles + + +def test_tiles_from_polygon(polygon): + """Get tiles from Polygon.""" + test_tiles = { + (9, 116, 544), + (9, 116, 545), + (9, 116, 546), + (9, 117, 540), + (9, 117, 541), + (9, 117, 542), + (9, 117, 543), + (9, 117, 544), + (9, 117, 545), + (9, 118, 536), + (9, 118, 537), + (9, 118, 538), + (9, 118, 539), + (9, 118, 540), + (9, 118, 541), + (9, 119, 535), + (9, 119, 536), + (9, 119, 537), + (9, 119, 538), + } + tp = TilePyramid("geodetic") + polygon_tiles = {tile.id for tile in tp.tiles_from_geom(polygon, 9)} + assert polygon_tiles == test_tiles + + +def test_tiles_from_multipolygon(multipolygon): + """Get tiles from MultiPolygon.""" + test_tiles = { + (9, 116, 544), + (9, 116, 545), + (9, 116, 546), + (9, 117, 540), + (9, 117, 541), + (9, 117, 542), + (9, 117, 543), + (9, 117, 544), + (9, 117, 545), + (9, 118, 536), + (9, 118, 537), + (9, 118, 538), + (9, 118, 539), + (9, 118, 540), + (9, 118, 541), + (9, 119, 535), + (9, 119, 536), + (9, 119, 537), + (9, 119, 538), + (9, 251, 604), + (9, 251, 605), + (9, 252, 604), + (9, 252, 605), + (9, 253, 605), + (9, 253, 606), + (9, 254, 605), + (9, 254, 606), + (9, 255, 606), + } + tp = TilePyramid("geodetic") + multipolygon_tiles = {tile.id for tile in tp.tiles_from_geom(multipolygon, 9)} + assert multipolygon_tiles == test_tiles + + +def test_tiles_from_point_batches(point): + """Get tile from point.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(point, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(point, zoom))) + + +def test_tiles_from_multipoint_batches(multipoint): + """Get tiles from multiple points.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(multipoint, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(multipoint, zoom))) + + +def test_tiles_from_linestring_batches(linestring): + """Get tiles from LineString.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(linestring, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(linestring, zoom))) + + +def test_tiles_from_multilinestring_batches(multilinestring): + """Get tiles from MultiLineString.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(multilinestring, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(multilinestring, zoom))) + + +def test_tiles_from_polygon_batches(polygon): + """Get tiles from Polygon.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(polygon, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(polygon, zoom))) + + +def test_tiles_from_multipolygon_batches(multipolygon): + """Get tiles from MultiPolygon.""" + tp = TilePyramid("geodetic") + zoom = 9 + tiles = 0 + gen = tp.tiles_from_geom(multipolygon, zoom, batch_by="row") + assert isinstance(gen, GeneratorType) + for row in gen: + assert isinstance(row, GeneratorType) + for tile in row: + tiles += 1 + assert isinstance(tile, Tile) + assert tiles + assert tiles == len(list(tp.tiles_from_geom(multipolygon, zoom))) diff --git a/test/test_tile.py b/test/test_tile.py new file mode 100644 index 0000000..52133df --- /dev/null +++ b/test/test_tile.py @@ -0,0 +1,200 @@ +"""Tile properties.""" + +import pytest +from affine import Affine + +from tilematrix import TilePyramid + + +def test_affine(): + """Affine output.""" + tp = TilePyramid("geodetic") + test_tiles = [(0, 0, 0), (1, 1, 1), (2, 2, 2)] + for tile_id in test_tiles: + tile = tp.tile(*tile_id) + test_affine = Affine( + tile.pixel_x_size, 0, tile.left, 0, -tile.pixel_y_size, tile.top + ) + assert tile.affine() == test_affine + # include pixelbuffer + pixelbuffer = 10 + test_tiles = [(1, 1, 1), (2, 2, 2), (3, 3, 3)] + for tile_id in test_tiles: + tile = tp.tile(*tile_id) + test_affine = Affine( + tile.pixel_x_size, + 0, + tile.bounds(pixelbuffer).left, + 0, + -tile.pixel_y_size, + tile.bounds(pixelbuffer).top, + ) + assert tile.affine(10) == test_affine + + +def test_get_parent(): + """Get parent Tile.""" + tp = TilePyramid("geodetic") + # default + tile = tp.tile(8, 100, 100) + assert tile.get_parent().id == (7, 50, 50) + # from top of pyramid + tile = tp.tile(0, 0, 0) + assert tile.get_parent() is None + + +def test_get_children(): + """Get Tile children.""" + # no metatiling + tp = TilePyramid("geodetic") + tile = tp.tile(8, 100, 100) + test_children = {(9, 200, 200), (9, 201, 200), (9, 200, 201), (9, 201, 201)} + children = {t.id for t in tile.get_children()} + assert test_children == children + + # 2 metatiling + tp = TilePyramid("geodetic", metatiling=2) + tile = tp.tile(0, 0, 0) + test_children = {(1, 0, 0), (1, 0, 1)} + children = {t.id for t in tile.get_children()} + assert test_children == children + + # 4 metatiling + tp = TilePyramid("geodetic", metatiling=4) + tile = tp.tile(0, 0, 0) + test_children = {(1, 0, 0)} + children = {t.id for t in tile.get_children()} + assert test_children == children + + +def test_get_neighbors(grid_definition_proj): + """Get Tile neighbors.""" + tp = TilePyramid("geodetic") + + # default + tile = tp.tile(8, 100, 100) + # 8 neighbors + test_neighbors = { + (8, 101, 100), + (8, 100, 101), + (8, 99, 100), + (8, 100, 99), + (8, 99, 101), + (8, 101, 101), + (8, 101, 99), + (8, 99, 99), + } + neighbors = {t.id for t in tile.get_neighbors()} + assert test_neighbors == neighbors + # 4 neighbors + test_neighbors = {(8, 101, 100), (8, 100, 101), (8, 99, 100), (8, 100, 99)} + neighbors = {t.id for t in tile.get_neighbors(connectedness=4)} + assert test_neighbors == neighbors + + # over antimeridian + tile = tp.tile(3, 1, 0) + # 8 neighbors + test_neighbors = { + (3, 0, 0), + (3, 1, 1), + (3, 2, 0), + (3, 1, 15), + (3, 0, 1), + (3, 2, 1), + (3, 2, 15), + (3, 0, 15), + } + neighbors = {t.id for t in tile.get_neighbors()} + assert test_neighbors == neighbors + # 4 neighbors + test_neighbors = {(3, 0, 0), (3, 1, 1), (3, 2, 0), (3, 1, 15)} + neighbors = {t.id for t in tile.get_neighbors(connectedness=4)} + assert test_neighbors == neighbors + + # tile has exactly two identical neighbors + tile = tp.tile(0, 0, 0) + test_tile = [(0, 0, 1)] + neighbors = [t.id for t in tile.get_neighbors()] + assert test_tile == neighbors + + # tile is alone at current zoom level + tp = TilePyramid("geodetic", metatiling=2) + tile = tp.tile(0, 0, 0) + neighbors = [t.id for t in tile.get_neighbors()] + assert neighbors == [] + + # wrong connectedness parameter + try: + tile.get_neighbors(connectedness="wrong_param") + raise Exception() + except ValueError: + pass + + # neighbors on non-global tilepyramids + tp = TilePyramid(grid_definition_proj) + zoom = 5 + max_col = tp.matrix_width(zoom) - 1 + tile = tp.tile(zoom=zoom, row=3, col=max_col) + # don't wrap around antimeridian, i.e. there are no tile neighbors + assert len(tile.get_neighbors()) == 5 + + +def test_intersecting(): + """Get intersecting Tiles from other TilePyramid.""" + tp_source = TilePyramid("geodetic", metatiling=2) + tp_target = TilePyramid("geodetic") + tile = tp_source.tile(5, 2, 2) + test_tiles = {(5, 4, 4), (5, 5, 4), (5, 4, 5), (5, 5, 5)} + intersecting_tiles = {t.id for t in tile.intersecting(tp_target)} + assert test_tiles == intersecting_tiles + + +def test_tile_compare(): + tp = TilePyramid("geodetic") + a = tp.tile(5, 5, 5) + b = tp.tile(5, 5, 5) + c = tp.tile(5, 5, 6) + assert a == b + assert a != c + assert b != c + assert a != "invalid type" + assert len(set([a, b, c])) == 2 + + +def test_tile_tuple(): + tp = TilePyramid("geodetic") + a = tp.tile(5, 5, 5) + assert tuple(a) == ( + 5, + 5, + 5, + ) + + +def test_deprecated(): + tp = TilePyramid("geodetic") + tile = tp.tile(5, 5, 5) + assert tile.srid + + +def test_invalid_id(): + tp = TilePyramid("geodetic") + # wrong id types + with pytest.raises(TypeError): + tp.tile(-1, 0, 0) + with pytest.raises(TypeError): + tp.tile(5, 0.7, 0) + with pytest.raises(TypeError): + tp.tile(10, 0, -11.56) + # column exceeds + with pytest.raises(ValueError): + tp.tile(5, 500, 0) + # row exceeds + with pytest.raises(ValueError): + tp.tile(5, 0, 500) + + +def test_tile_sizes(): + tp = TilePyramid("geodetic") + tile = tp.tile(5, 5, 5) + assert tile.x_size == tile.y_size diff --git a/test/test_tile_shapes.py b/test/test_tile_shapes.py new file mode 100644 index 0000000..de74dd3 --- /dev/null +++ b/test/test_tile_shapes.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python +"""Test Tile shapes.""" + +from tilematrix import TilePyramid + + +def test_simple_shapes(): + """Without metatiling & buffer.""" + # default + tile = TilePyramid("geodetic").tile(0, 0, 0) + assert tile.width == tile.height == 256 + assert tile.shape() == (256, 256) + # 512x512 + tile = TilePyramid("geodetic", tile_size=512).tile(0, 0, 0) + assert tile.width == tile.height == 512 + assert tile.shape() == (512, 512) + + +def test_geodetic_metatiling_shapes(): + """Metatile shapes.""" + # metatiling 2 + tp = TilePyramid("geodetic", metatiling=2) + assert tp.metatile_size == 512 + tile_shapes = { + (0, 0, 0): (256, 512), + (1, 0, 0): (512, 512), + (2, 0, 0): (512, 512), + (3, 0, 0): (512, 512), + (4, 0, 0): (512, 512), + (5, 0, 0): (512, 512), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.height == control_shape[0] + assert tile.width == control_shape[1] + assert tile.shape() == control_shape + + # metatiling 4 + tp = TilePyramid("geodetic", metatiling=4) + assert tp.metatile_size == 1024 + tile_shapes = { + (0, 0, 0): (256, 512), + (1, 0, 0): (512, 1024), + (2, 0, 0): (1024, 1024), + (3, 0, 0): (1024, 1024), + (4, 0, 0): (1024, 1024), + (5, 0, 0): (1024, 1024), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.height == control_shape[0] + assert tile.width == control_shape[1] + assert tile.shape() == control_shape + + # metatiling 8 + tp = TilePyramid("geodetic", metatiling=8) + assert tp.metatile_size == 2048 + tile_shapes = { + (0, 0, 0): (256, 512), + (1, 0, 0): (512, 1024), + (2, 0, 0): (1024, 2048), + (3, 0, 0): (2048, 2048), + (4, 0, 0): (2048, 2048), + (5, 0, 0): (2048, 2048), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.height == control_shape[0] + assert tile.width == control_shape[1] + assert tile.shape() == control_shape + + # metatiling 16 + tp = TilePyramid("geodetic", metatiling=16) + assert tp.metatile_size == 4096 + tile_shapes = { + (0, 0, 0): (256, 512), + (1, 0, 0): (512, 1024), + (2, 0, 0): (1024, 2048), + (3, 0, 0): (2048, 4096), + (4, 0, 0): (4096, 4096), + (5, 0, 0): (4096, 4096), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.height == control_shape[0] + assert tile.width == control_shape[1] + assert tile.shape() == control_shape + + +def test_geodetic_pixelbuffer_shapes(): + """Tile shapes using pixelbuffer.""" + tp = TilePyramid("geodetic") + pixelbuffer = 10 + tile_shapes = { + (0, 0, 0): (256, 276), # single tile at zoom 0 + (1, 0, 0): (266, 276), # top left + (2, 0, 0): (266, 276), # top left + (2, 0, 2): (266, 276), # top middle + (2, 0, 3): (266, 276), # top right + (2, 3, 0): (266, 276), # bottom left + (2, 3, 2): (266, 276), # bottom middle + (2, 3, 7): (266, 276), # bottom right + (3, 1, 0): (276, 276), # middle left + (3, 1, 1): (276, 276), # middle middle + (3, 1, 15): (276, 276), # middle right + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.shape(pixelbuffer) == control_shape + + +def test_geodetic_metatile_shapes(): + """Metatile shapes.""" + # metatiling 2 + tp = TilePyramid("geodetic", metatiling=2) + pixelbuffer = 10 + assert tp.metatile_size == 512 + tile_shapes = { + (0, 0, 0): (256, 532), + (1, 0, 0): (512, 532), + (2, 0, 0): (522, 532), + (3, 0, 0): (522, 532), + (4, 0, 0): (522, 532), + (5, 0, 0): (522, 532), + (5, 1, 1): (532, 532), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.shape(pixelbuffer) == control_shape + + # metatiling 4 + tp = TilePyramid("geodetic", metatiling=4) + assert tp.metatile_size == 1024 + tile_shapes = { + (0, 0, 0): (256, 532), + (1, 0, 0): (512, 1044), + (2, 0, 0): (1024, 1044), + (3, 0, 0): (1034, 1044), + (4, 0, 0): (1034, 1044), + (5, 0, 0): (1034, 1044), + (5, 1, 1): (1044, 1044), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.shape(pixelbuffer) == control_shape + + # metatiling 8 + tp = TilePyramid("geodetic", metatiling=8) + assert tp.metatile_size == 2048 + tile_shapes = { + (0, 0, 0): (256, 532), + (1, 0, 0): (512, 1044), + (2, 0, 0): (1024, 2068), + (3, 0, 0): (2048, 2068), + (4, 0, 0): (2058, 2068), + (5, 0, 0): (2058, 2068), + (5, 1, 1): (2068, 2068), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.shape(pixelbuffer) == control_shape + + # metatiling 16 + tp = TilePyramid("geodetic", metatiling=16) + assert tp.metatile_size == 4096 + tile_shapes = { + (0, 0, 0): (256, 532), + (1, 0, 0): (512, 1044), + (2, 0, 0): (1024, 2068), + (3, 0, 0): (2048, 4116), + (4, 0, 0): (4096, 4116), + (5, 0, 0): (4106, 4116), + (6, 1, 1): (4116, 4116), + } + for tile_id, control_shape in tile_shapes.items(): + tile = tp.tile(*tile_id) + assert tile.shape(pixelbuffer) == control_shape diff --git a/test/test_tilepyramid.py b/test/test_tilepyramid.py new file mode 100644 index 0000000..6d94cbf --- /dev/null +++ b/test/test_tilepyramid.py @@ -0,0 +1,376 @@ +"""TilePyramid creation.""" + +from types import GeneratorType + +import pytest +from shapely.geometry import Point, box +from shapely.ops import unary_union + +from tilematrix import TilePyramid, snap_bounds + + +def test_init(): + """Initialize TilePyramids.""" + for tptype in ["geodetic", "mercator"]: + assert TilePyramid(tptype) + with pytest.raises(ValueError): + TilePyramid("invalid") + with pytest.raises(ValueError): + TilePyramid() + assert hash(TilePyramid(tptype)) + + +def test_metatiling(): + """Metatiling setting.""" + for metatiling in [1, 2, 4, 8, 16]: + assert TilePyramid("geodetic", metatiling=metatiling) + try: + TilePyramid("geodetic", metatiling=5) + raise Exception() + except ValueError: + pass + + +def test_tile_size(): + """Tile sizes.""" + for tile_size in [128, 256, 512, 1024]: + tp = TilePyramid("geodetic", tile_size=tile_size) + assert tp.tile_size == tile_size + + +def test_intersect(): + """Get intersecting Tiles.""" + # same metatiling + tp = TilePyramid("geodetic") + intersect_tile = TilePyramid("geodetic").tile(5, 1, 1) + control = {(5, 1, 1)} + test_tiles = {tile.id for tile in tp.intersecting(intersect_tile)} + assert control == test_tiles + + # smaller metatiling + tp = TilePyramid("geodetic") + intersect_tile = TilePyramid("geodetic", metatiling=2).tile(5, 1, 1) + control = {(5, 2, 2), (5, 2, 3), (5, 3, 3), (5, 3, 2)} + test_tiles = {tile.id for tile in tp.intersecting(intersect_tile)} + assert control == test_tiles + + # bigger metatiling + tp = TilePyramid("geodetic", metatiling=2) + intersect_tile = TilePyramid("geodetic").tile(5, 1, 1) + control = {(5, 0, 0)} + test_tiles = {tile.id for tile in tp.intersecting(intersect_tile)} + assert control == test_tiles + intersect_tile = TilePyramid("geodetic").tile(4, 12, 31) + control = {(4, 6, 15)} + test_tiles = {tile.id for tile in tp.intersecting(intersect_tile)} + assert control == test_tiles + + # different CRSes + tp = TilePyramid("geodetic") + intersect_tile = TilePyramid("mercator").tile(5, 1, 1) + try: + test_tiles = {tile.id for tile in tp.intersecting(intersect_tile)} + raise Exception() + except ValueError: + pass + + +def test_tilepyramid_compare(grid_definition_proj, grid_definition_epsg): + """Comparison operators.""" + gproj, gepsg = grid_definition_proj, grid_definition_epsg + # predefined + assert TilePyramid("geodetic") == TilePyramid("geodetic") + assert TilePyramid("geodetic") != TilePyramid("geodetic", metatiling=2) + assert TilePyramid("geodetic") != TilePyramid("geodetic", tile_size=512) + assert TilePyramid("mercator") == TilePyramid("mercator") + assert TilePyramid("mercator") != TilePyramid("mercator", metatiling=2) + assert TilePyramid("mercator") != TilePyramid("mercator", tile_size=512) + # epsg based + assert TilePyramid(gepsg) == TilePyramid(gepsg) + assert TilePyramid(gepsg) != TilePyramid(gepsg, metatiling=2) + assert TilePyramid(gepsg) != TilePyramid(gepsg, tile_size=512) + # proj based + assert TilePyramid(gproj) == TilePyramid(gproj) + assert TilePyramid(gproj) != TilePyramid(gproj, metatiling=2) + assert TilePyramid(gproj) != TilePyramid(gproj, tile_size=512) + # altered bounds + abounds = dict(**gproj) + abounds.update(bounds=(-5000000.0, -5000000.0, 5000000.0, 5000000.0)) + assert TilePyramid(abounds) == TilePyramid(abounds) + assert TilePyramid(gproj) != TilePyramid(abounds) + # other type + assert TilePyramid("geodetic") != "string" + + +def test_grid_compare(grid_definition_proj, grid_definition_epsg): + """Comparison operators.""" + gproj, gepsg = grid_definition_proj, grid_definition_epsg + # predefined + assert TilePyramid("geodetic").grid == TilePyramid("geodetic").grid + assert TilePyramid("geodetic").grid == TilePyramid("geodetic", metatiling=2).grid + assert TilePyramid("geodetic").grid == TilePyramid("geodetic", tile_size=512).grid + assert TilePyramid("mercator").grid == TilePyramid("mercator").grid + assert TilePyramid("mercator").grid == TilePyramid("mercator", metatiling=2).grid + assert TilePyramid("mercator").grid == TilePyramid("mercator", tile_size=512).grid + # epsg based + assert TilePyramid(gepsg).grid == TilePyramid(gepsg).grid + assert TilePyramid(gepsg).grid == TilePyramid(gepsg, metatiling=2).grid + assert TilePyramid(gepsg).grid == TilePyramid(gepsg, tile_size=512).grid + # proj based + assert TilePyramid(gproj).grid == TilePyramid(gproj).grid + assert TilePyramid(gproj).grid == TilePyramid(gproj, metatiling=2).grid + assert TilePyramid(gproj).grid == TilePyramid(gproj, tile_size=512).grid + # altered bounds + abounds = dict(**gproj) + abounds.update(bounds=(-5000000.0, -5000000.0, 5000000.0, 5000000.0)) + assert TilePyramid(abounds).grid == TilePyramid(abounds).grid + assert TilePyramid(gproj).grid != TilePyramid(abounds).grid + + +def test_tile_from_xy(): + tp = TilePyramid("geodetic") + zoom = 5 + + # point inside tile + p_in = (0.5, 0.5, zoom) + control_in = [ + ((5, 15, 32), "rb"), + ((5, 15, 32), "lb"), + ((5, 15, 32), "rt"), + ((5, 15, 32), "lt"), + ] + for tile_id, on_edge_use in control_in: + tile = tp.tile_from_xy(*p_in, on_edge_use=on_edge_use) + assert tile.id == tile_id + assert Point(p_in[0], p_in[1]).within(tile.bbox()) + + # point is on tile edge + p_edge = (0, 0, zoom) + control_edge = [ + ((5, 16, 32), "rb"), + ((5, 16, 31), "lb"), + ((5, 15, 32), "rt"), + ((5, 15, 31), "lt"), + ] + for tile_id, on_edge_use in control_edge: + tile = tp.tile_from_xy(*p_edge, on_edge_use=on_edge_use) + assert tile.id == tile_id + assert Point(p_edge[0], p_edge[1]).touches(tile.bbox()) + + with pytest.raises(ValueError): + tp.tile_from_xy(180, -90, zoom, on_edge_use="rb") + with pytest.raises(ValueError): + tp.tile_from_xy(180, -90, zoom, on_edge_use="lb") + + tile = tp.tile_from_xy(180, -90, zoom, on_edge_use="rt") + assert tile.id == (5, 31, 0) + tile = tp.tile_from_xy(180, -90, zoom, on_edge_use="lt") + assert tile.id == (5, 31, 63) + + with pytest.raises(TypeError): + tp.tile_from_xy(-180, 90, zoom, on_edge_use="lt") + with pytest.raises(TypeError): + tp.tile_from_xy(-180, 90, zoom, on_edge_use="rt") + + tile = tp.tile_from_xy(-180, 90, zoom, on_edge_use="rb") + assert tile.id == (5, 0, 0) + tile = tp.tile_from_xy(-180, 90, zoom, on_edge_use="lb") + assert tile.id == (5, 0, 63) + + with pytest.raises(ValueError): + tp.tile_from_xy(-180, 90, zoom, on_edge_use="invalid") + + +def test_tiles_from_bounds(grid_definition_proj): + # global pyramids + tp = TilePyramid("geodetic") + parent = tp.tile(8, 5, 5) + from_bounds = set([t.id for t in tp.tiles_from_bounds(parent.bounds(), 9)]) + children = set([t.id for t in parent.get_children()]) + assert from_bounds == children + # non-global pyramids + tp = TilePyramid(grid_definition_proj) + parent = tp.tile(8, 0, 0) + from_bounds = set([t.id for t in tp.tiles_from_bounds(parent.bounds(), 9)]) + children = set([t.id for t in parent.get_children()]) + assert from_bounds == children + + +def test_tiles_from_bounds_batch_by_row(): + tp = TilePyramid("geodetic") + bounds = (0, 0, 90, 90) + zoom = 8 + + tiles = tp.tiles_from_bounds(bounds, zoom, batch_by="row") + assert isinstance(tiles, GeneratorType) + assert list(tiles) + + previous_row = None + tiles = 0 + for tile_row in tp.tiles_from_bounds(bounds, zoom, batch_by="row"): + assert isinstance(tile_row, GeneratorType) + previous_tile = None + for tile in tile_row: + tiles += 1 + if previous_row is None: + if previous_tile is not None: + assert tile.col == previous_tile.col + 1 + else: + if previous_tile is not None: + assert tile.col == previous_tile.col + 1 + assert tile.row == previous_tile.row + assert tile.row == previous_row + 1 + + previous_tile = tile + + previous_row = tile.row + + assert tiles == len(list(tp.tiles_from_bounds(bounds, zoom))) + + +def test_tiles_from_bounds_batch_by_column(): + tp = TilePyramid("geodetic") + bounds = (0, 0, 90, 90) + zoom = 8 + + tiles = tp.tiles_from_bounds(bounds, zoom, batch_by="column") + assert isinstance(tiles, GeneratorType) + assert list(tiles) + + previous_column = None + tiles = 0 + for tile_column in tp.tiles_from_bounds(bounds, zoom, batch_by="column"): + assert isinstance(tile_column, GeneratorType) + previous_tile = None + for tile in tile_column: + tiles += 1 + if previous_column is None: + if previous_tile is not None: + assert tile.row == previous_tile.row + 1 + else: + if previous_tile is not None: + assert tile.row == previous_tile.row + 1 + assert tile.col == previous_tile.col + assert tile.col == previous_column + 1 + + previous_tile = tile + + previous_column = tile.col + + assert tiles == len(list(tp.tiles_from_bounds(bounds, zoom))) + + +def test_tiles_from_bounds_batch_by_row_antimeridian_bounds(): + tp = TilePyramid("geodetic") + bounds = (0, 0, 185, 95) + zoom = 8 + + tiles = tp.tiles_from_bounds(bounds, zoom, batch_by="row") + assert isinstance(tiles, GeneratorType) + assert list(tiles) + + previous_row = None + tiles = 0 + for tile_row in tp.tiles_from_bounds(bounds, zoom, batch_by="row"): + assert isinstance(tile_row, GeneratorType) + previous_tile = None + for tile in tile_row: + tiles += 1 + if previous_row is None: + if previous_tile is not None: + assert tile.col > previous_tile.col + else: + if previous_tile is not None: + assert tile.col > previous_tile.col + assert tile.row == previous_tile.row + assert tile.row > previous_row + + previous_tile = tile + + previous_row = tile.row + + assert tiles == len(list(tp.tiles_from_bounds(bounds, zoom))) + + +def test_tiles_from_bounds_batch_by_row_both_antimeridian_bounds(): + tp = TilePyramid("geodetic") + bounds = (-185, 0, 185, 95) + zoom = 8 + + tiles = tp.tiles_from_bounds(bounds, zoom, batch_by="row") + assert isinstance(tiles, GeneratorType) + assert list(tiles) + + previous_row = None + tiles = 0 + for tile_row in tp.tiles_from_bounds(bounds, zoom, batch_by="row"): + assert isinstance(tile_row, GeneratorType) + previous_tile = None + for tile in tile_row: + tiles += 1 + if previous_row is None: + if previous_tile is not None: + assert tile.col == previous_tile.col + 1 + else: + if previous_tile is not None: + assert tile.col == previous_tile.col + 1 + assert tile.row == previous_tile.row + assert tile.row == previous_row + 1 + + previous_tile = tile + + previous_row = tile.row + + assert tiles == len(list(tp.tiles_from_bounds(bounds, zoom))) + + +def test_tiles_from_geom_exact(tile_bounds_polygon): + tp = TilePyramid("geodetic") + zoom = 3 + + tiles = len(list(tp.tiles_from_geom(tile_bounds_polygon, zoom))) + assert tiles == 4 + tiles = 0 + for batch in tp.tiles_from_geom(tile_bounds_polygon, zoom, batch_by="row"): + tiles += len(list(batch)) + assert tiles == 4 + + exact_tiles = len(list(tp.tiles_from_geom(tile_bounds_polygon, zoom, exact=True))) + assert exact_tiles == 3 + tiles = 0 + for batch in tp.tiles_from_geom( + tile_bounds_polygon, zoom, batch_by="row", exact=True + ): + tiles += len(list(batch)) + assert tiles == 3 + + +def test_snap_bounds(): + bounds = (0, 1, 2, 3) + tp = TilePyramid("geodetic") + zoom = 8 + + snapped = snap_bounds(bounds=bounds, tile_pyramid=tp, zoom=zoom) + control = unary_union( + [tile.bbox() for tile in tp.tiles_from_bounds(bounds, zoom)] + ).bounds + assert snapped == control + + pixelbuffer = 10 + snapped = snap_bounds( + bounds=bounds, tile_pyramid=tp, zoom=zoom, pixelbuffer=pixelbuffer + ) + control = unary_union( + [tile.bbox(pixelbuffer) for tile in tp.tiles_from_bounds(bounds, zoom)] + ).bounds + assert snapped == control + + +def test_deprecated(): + tp = TilePyramid("geodetic") + assert tp.type + assert tp.srid + assert tp.tile_x_size(0) + assert tp.tile_y_size(0) + assert tp.tile_height(0) + assert tp.tile_width(0) diff --git a/tilematrix/__init__.py b/tilematrix/__init__.py new file mode 100644 index 0000000..8e16239 --- /dev/null +++ b/tilematrix/__init__.py @@ -0,0 +1,24 @@ +"""Main package entry point.""" + +from ._conf import PYRAMID_PARAMS +from ._funcs import clip_geometry_to_srs_bounds, snap_bounds, validate_zoom +from ._grid import GridDefinition +from ._tile import Tile +from ._tilepyramid import TilePyramid +from ._types import Bounds, Shape, TileIndex + +__all__ = [ + "Bounds", + "clip_geometry_to_srs_bounds", + "GridDefinition", + "PYRAMID_PARAMS", + "Shape", + "snap_bounds", + "TilePyramid", + "Tile", + "TileIndex", + "validate_zoom", +] + + +__version__ = "2023.12.0" diff --git a/tilematrix/_conf.py b/tilematrix/_conf.py new file mode 100644 index 0000000..4eacdac --- /dev/null +++ b/tilematrix/_conf.py @@ -0,0 +1,28 @@ +"""Package configuration parameters.""" + +# round coordinates +ROUND = 20 + +# bounds ratio vs shape ratio uncertainty +DELTA = 1e-6 + +# supported pyramid types +PYRAMID_PARAMS = { + "geodetic": { + "shape": (1, 2), # tile rows and columns at zoom level 0 + "bounds": (-180.0, -90.0, 180.0, 90.0), # pyramid bounds in pyramid CRS + "is_global": True, # if false, no antimeridian handling + "srs": {"epsg": 4326}, # EPSG code for CRS + }, + "mercator": { + "shape": (1, 1), + "bounds": ( + -20037508.3427892, + -20037508.3427892, + 20037508.3427892, + 20037508.3427892, + ), + "is_global": True, + "srs": {"epsg": 3857}, + }, +} diff --git a/tilematrix/_funcs.py b/tilematrix/_funcs.py new file mode 100644 index 0000000..ee5eceb --- /dev/null +++ b/tilematrix/_funcs.py @@ -0,0 +1,263 @@ +"""Helper functions.""" + +from itertools import product + +from rasterio.crs import CRS +from shapely.affinity import translate +from shapely.geometry import GeometryCollection, Polygon, box +from shapely.ops import unary_union +from shapely.prepared import prep + +from ._conf import DELTA, ROUND +from ._types import Bounds, Shape + + +def validate_zoom(zoom): + if not isinstance(zoom, int): + raise TypeError("zoom must be an integer") + if zoom < 0: + raise ValueError("zoom must be greater or equal 0") + + +def clip_geometry_to_srs_bounds(geometry, pyramid, multipart=False): + """ + Clip input geometry to SRS bounds of given TilePyramid. + + If geometry passes the antimeridian, it will be split up in a multipart + geometry and shifted to within the SRS boundaries. + Note: geometry SRS must be the TilePyramid SRS! + + - geometry: any shapely geometry + - pyramid: a TilePyramid object + - multipart: return list of geometries instead of a GeometryCollection + """ + if not geometry.is_valid: + raise ValueError("invalid geometry given") + pyramid_bbox = box(*pyramid.bounds) + + # Special case for global tile pyramids if geometry extends over tile + # pyramid boundaries (such as the antimeridian). + if pyramid.is_global and not geometry.within(pyramid_bbox): + inside_geom = geometry.intersection(pyramid_bbox) + outside_geom = geometry.difference(pyramid_bbox) + # shift outside geometry so it lies within SRS bounds + if hasattr(outside_geom, "geoms"): + outside_geoms = outside_geom.geoms + else: + outside_geoms = [outside_geom] + all_geoms = [inside_geom] + for geom in outside_geoms: + geom_bounds = Bounds(*geom.bounds) + if geom_bounds.left < pyramid.left: + geom = translate(geom, xoff=2 * pyramid.right) + elif geom_bounds.right > pyramid.right: + geom = translate(geom, xoff=-2 * pyramid.right) + all_geoms.append(geom) + if multipart: + return all_geoms + else: + return GeometryCollection(all_geoms) + + else: + if multipart: + return [geometry] + else: + return geometry + + +def snap_bounds(bounds=None, tile_pyramid=None, zoom=None, pixelbuffer=0): + """ + Extend bounds to be aligned with union of tile bboxes. + + - bounds: (left, bottom, right, top) + - tile_pyramid: a TilePyramid object + - zoom: target zoom level + - pixelbuffer: apply pixelbuffer + """ + bounds = Bounds(*bounds) + validate_zoom(zoom) + lb = _tile_from_xy(tile_pyramid, bounds.left, bounds.bottom, zoom, on_edge_use="rt") + rt = _tile_from_xy(tile_pyramid, bounds.right, bounds.top, zoom, on_edge_use="lb") + left, bottom, _, _ = lb.bounds(pixelbuffer) + _, _, right, top = rt.bounds(pixelbuffer) + return Bounds(left, bottom, right, top) + + +def _verify_shape_bounds(shape, bounds): + """Verify that shape corresponds to bounds apect ratio.""" + if not isinstance(shape, (tuple, list)) or len(shape) != 2: + raise TypeError( + "shape must be a tuple or list with two elements: %s" % str(shape) + ) + if not isinstance(bounds, (tuple, list)) or len(bounds) != 4: + raise TypeError( + "bounds must be a tuple or list with four elements: %s" % str(bounds) + ) + shape = Shape(*shape) + bounds = Bounds(*bounds) + shape_ratio = shape.width / shape.height + bounds_ratio = (bounds.right - bounds.left) / (bounds.top - bounds.bottom) + if abs(shape_ratio - bounds_ratio) > DELTA: + min_length = min( + [ + (bounds.right - bounds.left) / shape.width, + (bounds.top - bounds.bottom) / shape.height, + ] + ) + proposed_bounds = Bounds( + bounds.left, + bounds.bottom, + bounds.left + shape.width * min_length, + bounds.bottom + shape.height * min_length, + ) + raise ValueError( + "shape ratio (%s) must equal bounds ratio (%s); try %s" + % (shape_ratio, bounds_ratio, proposed_bounds) + ) + + +def _get_crs(srs): + if not isinstance(srs, dict): + raise TypeError("'srs' must be a dictionary") + if "wkt" in srs: + return CRS().from_wkt(srs["wkt"]) + elif "epsg" in srs: + return CRS().from_epsg(srs["epsg"]) + elif "proj" in srs: + return CRS().from_string(srs["proj"]) + else: + raise TypeError("provide either 'wkt', 'epsg' or 'proj' definition") + + +def _tile_intersecting_tilepyramid(tile, tp): + """Return all tiles from tilepyramid intersecting with tile.""" + if tile.tp.grid != tp.grid: + raise ValueError("Tile and TilePyramid source grids must be the same.") + tile_metatiling = tile.tile_pyramid.metatiling + pyramid_metatiling = tp.metatiling + multiplier = tile_metatiling / pyramid_metatiling + if tile_metatiling > pyramid_metatiling: + return [ + tp.tile( + tile.zoom, + int(multiplier) * tile.row + row_offset, + int(multiplier) * tile.col + col_offset, + ) + for row_offset, col_offset in product( + range(int(multiplier)), range(int(multiplier)) + ) + ] + elif tile_metatiling < pyramid_metatiling: + return [ + tp.tile(tile.zoom, int(multiplier * tile.row), int(multiplier * tile.col)) + ] + else: + return [tp.tile(*tile.id)] + + +def _global_tiles_from_bounds(tp, bounds, zoom, batch_by=None): + """Return also Tiles if bounds cross the antimeridian.""" + + # clip to tilepyramid top and bottom bounds + left, right = bounds.left, bounds.right + top = min([tp.top, bounds.top]) + bottom = max([tp.bottom, bounds.bottom]) + + # special case if bounds cross the antimeridian + if left < tp.left or right > tp.right: + # shift overlap to valid side of antimeridian + # create MultiPolygon + # yield tiles by + bounds_geoms = [] + + # tiles west of antimeridian + if left < tp.left: + # move outside part of bounding box to the valid side of antimeridian + bounds_geoms.append(box(left + (tp.right - tp.left), bottom, tp.right, top)) + bounds_geoms.append( + # remaining part of bounding box + box(tp.left, bottom, min([right, tp.right]), top) + ) + + # tiles east of antimeridian + if right > tp.right: + # move outside part of bounding box to the valid side of antimeridian + bounds_geoms.append(box(tp.left, bottom, right - (tp.right - tp.left), top)) + bounds_geoms.append(box(max([left, tp.left]), bottom, tp.right, top)) + + bounds_geom = unary_union(bounds_geoms).buffer(0) + bounds_geom_prep = prep(bounds_geom) + + # if union of bounding boxes is a multipart geometry, do some costly checks to be able + # to yield in batches + if bounds_geom.geom_type.lower().startswith("multi"): + if batch_by: + for batch in _tiles_from_cleaned_bounds( + tp, bounds_geom.bounds, zoom, batch_by + ): + yield ( + tile + for tile in batch + if bounds_geom_prep.intersects(tile.bbox()) + ) + else: + for tile in _tiles_from_cleaned_bounds(tp, bounds_geom.bounds, zoom): + if bounds_geom_prep.intersects(tile.bbox()): + yield tile + return + + # else, continue with cleaned bounds + bounds = bounds_geom.bounds + + # yield using cleaned bounds + yield from _tiles_from_cleaned_bounds(tp, bounds, zoom, batch_by=batch_by) + + +def _tiles_from_cleaned_bounds(tp, bounds, zoom, batch_by=None): + """Return all tiles intersecting with bounds.""" + bounds = Bounds(*bounds) + lb = _tile_from_xy(tp, bounds.left, bounds.bottom, zoom, on_edge_use="rt") + rt = _tile_from_xy(tp, bounds.right, bounds.top, zoom, on_edge_use="lb") + row_range = range(rt.row, lb.row + 1) + col_range = range(lb.col, rt.col + 1) + if batch_by is None: + for row, col in product(row_range, col_range): + yield tp.tile(zoom, row, col) + elif batch_by == "row": + for row in row_range: + yield (tp.tile(zoom, row, col) for col in col_range) + elif batch_by == "column": + for col in col_range: + yield (tp.tile(zoom, row, col) for row in row_range) + else: # pragma: no cover + raise ValueError("'batch_by' must either be None, 'row' or 'column'.") + + +def _tile_from_xy(tp, x, y, zoom, on_edge_use="rb"): + # determine row + tile_y_size = round(tp.pixel_y_size(zoom) * tp.tile_size * tp.metatiling, ROUND) + row = int((tp.top - y) / tile_y_size) + if on_edge_use in ["rt", "lt"] and (tp.top - y) % tile_y_size == 0.0: + row -= 1 + + # determine column + tile_x_size = round(tp.pixel_x_size(zoom) * tp.tile_size * tp.metatiling, ROUND) + col = int((x - tp.left) / tile_x_size) + if on_edge_use in ["lb", "lt"] and (x - tp.left) % tile_x_size == 0.0: + col -= 1 + + # handle Antimeridian wrapping + if tp.is_global: + # left side + if col == -1: + col = tp.matrix_width(zoom) - 1 + # right side + elif col >= tp.matrix_width(zoom): + col = col % tp.matrix_width(zoom) + + try: + return tp.tile(zoom, row, col) + except ValueError as e: + raise ValueError( + "on_edge_use '%s' results in an invalid tile: %s" % (on_edge_use, e) + ) diff --git a/tilematrix/_grid.py b/tilematrix/_grid.py new file mode 100644 index 0000000..03376e3 --- /dev/null +++ b/tilematrix/_grid.py @@ -0,0 +1,107 @@ +import warnings + +from ._conf import PYRAMID_PARAMS +from ._funcs import _get_crs, _verify_shape_bounds +from ._types import Bounds, Shape + + +class GridDefinition(object): + """Object representing the tile pyramid source grid.""" + + def __init__( + self, grid=None, shape=None, bounds=None, srs=None, is_global=False, **kwargs + ): + if isinstance(grid, str) and grid in PYRAMID_PARAMS: + self.type = grid + self.shape = Shape(*PYRAMID_PARAMS[grid]["shape"]) + self.bounds = Bounds(*PYRAMID_PARAMS[grid]["bounds"]) + self.is_global = PYRAMID_PARAMS[grid]["is_global"] + self.crs = _get_crs(PYRAMID_PARAMS[grid]["srs"]) + self.left, self.bottom, self.right, self.top = self.bounds + elif grid is None or grid == "custom": + for i in ["proj", "epsg"]: + if i in kwargs: + srs = {i: kwargs[i]} if srs is None else srs + warnings.warn( + DeprecationWarning( + "'%s' should be packed into a dictionary and passed to " + "'srs'" % i + ) + ) + self.type = "custom" + _verify_shape_bounds(shape=shape, bounds=bounds) + self.shape = Shape(*shape) + self.bounds = Bounds(*bounds) + self.is_global = is_global + self.crs = _get_crs(srs) + self.left, self.bottom, self.right, self.top = self.bounds + # check if parameters match with default grid type + for default_grid_name in PYRAMID_PARAMS: + default_grid = GridDefinition(default_grid_name) + if self.__eq__(default_grid): + self.type = default_grid_name + elif isinstance(grid, dict): + if "type" in grid: # pragma: no cover + warnings.warn( + DeprecationWarning("'type' is deprecated and should be 'grid'") + ) + if "grid" not in grid: + grid["grid"] = grid.pop("type") + self.__init__(**grid) + elif isinstance(grid, GridDefinition): + self.__init__(**grid.to_dict()) + else: + raise ValueError("invalid grid definition: %s" % grid) + + @property + def srid(self): + warnings.warn(DeprecationWarning("'srid' attribute is deprecated")) + return self.crs.to_epsg() + + def to_dict(self): + return dict( + bounds=self.bounds, + is_global=self.is_global, + shape=self.shape, + srs=dict(wkt=self.crs.to_wkt()), + grid=self.type, + ) + + def from_dict(config_dict): + return GridDefinition(**config_dict) + + def __eq__(self, other): + return ( + isinstance(other, self.__class__) + and self.shape == other.shape + and self.bounds == other.bounds + and self.is_global == other.is_global + and self.crs == other.crs + ) + + def __ne__(self, other): + return not self.__eq__(other) + + def __repr__(self): + if self.type in PYRAMID_PARAMS: + return 'GridDefinition("%s")' % self.type + else: + return ( + "GridDefinition(" + '"%s", ' + "shape=%s, " + "bounds=%s, " + "is_global=%s, " + "srs=%s" + ")" + % ( + self.type, + tuple(self.shape), + tuple(self.bounds), + self.is_global, + self.crs, + ) + ) + + def __hash__(self): + return hash(repr(self)) diff --git a/tilematrix/_tile.py b/tilematrix/_tile.py new file mode 100644 index 0000000..66ab649 --- /dev/null +++ b/tilematrix/_tile.py @@ -0,0 +1,304 @@ +"""Tile class.""" + +import warnings + +from affine import Affine +from shapely.geometry import box + +from ._conf import ROUND +from ._funcs import _tile_intersecting_tilepyramid +from ._types import Bounds, Shape, TileIndex + + +class Tile(object): + """ + A Tile is a square somewhere on Earth. + + Each Tile can be identified with the zoom, row, column index in a + TilePyramid. + + Some tile functions can accept a tile buffer in pixels (pixelbuffer). A + pixelbuffer value of e.g. 1 will extend the tile boundaries by 1 pixel. + """ + + def __init__(self, tile_pyramid, zoom, row, col): + """Initialize Tile.""" + self.tile_pyramid = tile_pyramid + self.tp = tile_pyramid + self.crs = tile_pyramid.crs + self.zoom = zoom + self.row = row + self.col = col + self.is_valid() + self.index = self.id = TileIndex(zoom, row, col) + self.pixel_x_size = self.tile_pyramid.pixel_x_size(self.zoom) + self.pixel_y_size = self.tile_pyramid.pixel_y_size(self.zoom) + # base SRID size without pixelbuffer + self._base_srid_size = Shape( + height=self.pixel_y_size * self.tp.tile_size * self.tp.metatiling, + width=self.pixel_x_size * self.tp.tile_size * self.tp.metatiling, + ) + # base bounds not accounting for pixelbuffers but metatiles are clipped to + # TilePyramid bounds + self._top = round(self.tp.top - (self.row * self._base_srid_size.height), ROUND) + self._bottom = max([self._top - self._base_srid_size.height, self.tp.bottom]) + self._left = round( + self.tp.left + (self.col * self._base_srid_size.width), ROUND + ) + self._right = min([self._left + self._base_srid_size.width, self.tp.right]) + # base shape without pixelbuffer + self._base_shape = Shape( + height=int(round((self._top - self._bottom) / self.pixel_y_size, 0)), + width=int(round((self._right - self._left) / self.pixel_x_size, 0)), + ) + + @property + def left(self): + return self.bounds().left + + @property + def bottom(self): + return self.bounds().bottom + + @property + def right(self): + return self.bounds().right + + @property + def top(self): + return self.bounds().top + + @property + def srid(self): + warnings.warn(DeprecationWarning("'srid' attribute is deprecated")) + return self.tp.grid.srid + + @property + def width(self): + """Calculate Tile width in pixels.""" + return self.shape().width + + @property + def height(self): + """Calculate Tile height in pixels.""" + return self.shape().height + + @property + def x_size(self): + """Width of tile in SRID units at zoom level.""" + return self.right - self.left + + @property + def y_size(self): + """Height of tile in SRID units at zoom level.""" + return self.top - self.bottom + + def bounds(self, pixelbuffer=0): + """ + Return Tile boundaries. + + - pixelbuffer: tile buffer in pixels + """ + left = self._left + bottom = self._bottom + right = self._right + top = self._top + if pixelbuffer: + offset = self.pixel_x_size * float(pixelbuffer) + left -= offset + bottom -= offset + right += offset + top += offset + # on global grids clip at northern and southern TilePyramid bound + if self.tp.grid.is_global: + top = min([top, self.tile_pyramid.top]) + bottom = max([bottom, self.tile_pyramid.bottom]) + return Bounds(left, bottom, right, top) + + def bbox(self, pixelbuffer=0): + """ + Return Tile bounding box. + + - pixelbuffer: tile buffer in pixels + """ + return box(*self.bounds(pixelbuffer=pixelbuffer)) + + def affine(self, pixelbuffer=0): + """ + Return an Affine object of tile. + + - pixelbuffer: tile buffer in pixels + """ + return Affine( + self.pixel_x_size, + 0, + self.bounds(pixelbuffer).left, + 0, + -self.pixel_y_size, + self.bounds(pixelbuffer).top, + ) + + def shape(self, pixelbuffer=0): + """ + Return a tuple of tile height and width. + + - pixelbuffer: tile buffer in pixels + """ + # apply pixelbuffers + height = self._base_shape.height + 2 * pixelbuffer + width = self._base_shape.width + 2 * pixelbuffer + if pixelbuffer and self.tp.grid.is_global: + # on first and last row, remove pixelbuffer on top or bottom + matrix_height = self.tile_pyramid.matrix_height(self.zoom) + if matrix_height == 1: + height = self._base_shape.height + elif self.row in [0, matrix_height - 1]: + height = self._base_shape.height + pixelbuffer + return Shape(height=height, width=width) + + def is_valid(self): + """Return True if tile is available in tile pyramid.""" + if not all( + [ + isinstance(self.zoom, int), + self.zoom >= 0, + isinstance(self.row, int), + self.row >= 0, + isinstance(self.col, int), + self.col >= 0, + ] + ): + raise TypeError("zoom, col and row must be integers >= 0") + cols = self.tile_pyramid.matrix_width(self.zoom) + rows = self.tile_pyramid.matrix_height(self.zoom) + if self.col >= cols: + raise ValueError("col (%s) exceeds matrix width (%s)" % (self.col, cols)) + if self.row >= rows: + raise ValueError("row (%s) exceeds matrix height (%s)" % (self.row, rows)) + return True + + def get_parent(self): + """Return tile from previous zoom level.""" + return ( + None + if self.zoom == 0 + else self.tile_pyramid.tile(self.zoom - 1, self.row // 2, self.col // 2) + ) + + def get_children(self): + """Return tiles from next zoom level.""" + next_zoom = self.zoom + 1 + return [ + self.tile_pyramid.tile( + next_zoom, self.row * 2 + row_offset, self.col * 2 + col_offset + ) + for row_offset, col_offset in [ + (0, 0), # top left + (0, 1), # top right + (1, 1), # bottom right + (1, 0), # bottom left + ] + if all( + [ + self.row * 2 + row_offset < self.tp.matrix_height(next_zoom), + self.col * 2 + col_offset < self.tp.matrix_width(next_zoom), + ] + ) + ] + + def get_neighbors(self, connectedness=8): + """ + Return tile neighbors. + + Tile neighbors are unique, i.e. in some edge cases, where both the left + and right neighbor wrapped around the antimeridian is the same. Also, + neighbors ouside the northern and southern TilePyramid boundaries are + excluded, because they are invalid. + + ------------- + | 8 | 1 | 5 | + ------------- + | 4 | x | 2 | + ------------- + | 7 | 3 | 6 | + ------------- + + - connectedness: [4 or 8] return four direct neighbors or all eight. + """ + if connectedness not in [4, 8]: + raise ValueError("only connectedness values 8 or 4 are allowed") + + unique_neighbors = {} + # 4-connected neighborsfor pyramid + matrix_offsets = [ + (-1, 0), # 1: above + (0, 1), # 2: right + (1, 0), # 3: below + (0, -1), # 4: left + ] + if connectedness == 8: + matrix_offsets.extend( + [ + (-1, 1), # 5: above right + (1, 1), # 6: below right + (1, -1), # 7: below left + (-1, -1), # 8: above left + ] + ) + + for row_offset, col_offset in matrix_offsets: + new_row = self.row + row_offset + new_col = self.col + col_offset + # omit if row is outside of tile matrix + if new_row < 0 or new_row >= self.tp.matrix_height(self.zoom): + continue + # wrap around antimeridian if new column is outside of tile matrix + if new_col < 0: + if not self.tp.is_global: + continue + new_col = self.tp.matrix_width(self.zoom) + new_col + elif new_col >= self.tp.matrix_width(self.zoom): + if not self.tp.is_global: + continue + new_col -= self.tp.matrix_width(self.zoom) + # omit if new tile is current tile + if new_row == self.row and new_col == self.col: + continue + # create new tile + unique_neighbors[(new_row, new_col)] = self.tp.tile( + self.zoom, new_row, new_col + ) + + return unique_neighbors.values() + + def intersecting(self, tilepyramid): + """ + Return all Tiles from intersecting TilePyramid. + + This helps translating between TilePyramids with different metatiling + settings. + + - tilepyramid: a TilePyramid object + """ + return _tile_intersecting_tilepyramid(self, tilepyramid) + + def __eq__(self, other): + return ( + isinstance(other, self.__class__) + and self.tp == other.tp + and self.id == other.id + ) + + def __ne__(self, other): + return not self.__eq__(other) + + def __repr__(self): + return "Tile(%s, %s)" % (self.id, self.tp) + + def __hash__(self): + return hash(repr(self)) + + def __iter__(self): + yield self.zoom + yield self.row + yield self.col diff --git a/tilematrix/_tilepyramid.py b/tilematrix/_tilepyramid.py new file mode 100644 index 0000000..d85c11c --- /dev/null +++ b/tilematrix/_tilepyramid.py @@ -0,0 +1,326 @@ +"""Handling tile pyramids.""" + +import math +import warnings + +from shapely.prepared import prep + +from ._conf import ROUND +from ._funcs import ( + _global_tiles_from_bounds, + _tile_from_xy, + _tile_intersecting_tilepyramid, + _tiles_from_cleaned_bounds, + clip_geometry_to_srs_bounds, + validate_zoom, +) +from ._grid import GridDefinition +from ._tile import Tile +from ._types import Bounds + + +class TilePyramid(object): + """ + A Tile pyramid is a collection of tile matrices for different zoom levels. + + TilePyramids can be defined either for the Web Mercator or the geodetic + projection. A TilePyramid holds tiles for different discrete zoom levels. + Each zoom level is a 2D tile matrix, thus every Tile can be defined by + providing the zoom level as well as the row and column of the tile matrix. + + - projection: one of "geodetic" or "mercator" + - tile_size: target pixel size of each tile + - metatiling: tile size mulipilcation factor, must be one of 1, 2, 4, 8 or + 16. + """ + + def __init__(self, grid=None, tile_size=256, metatiling=1): + """Initialize TilePyramid.""" + if grid is None: + raise ValueError("grid definition required") + _metatiling_opts = [2**x for x in range(10)] + if metatiling not in _metatiling_opts: + raise ValueError(f"metatling must be one of {_metatiling_opts}") + # get source grid parameters + self.grid = GridDefinition(grid) + self.bounds = self.grid.bounds + self.left, self.bottom, self.right, self.top = self.bounds + self.crs = self.grid.crs + self.is_global = self.grid.is_global + self.metatiling = metatiling + # size in pixels + self.tile_size = tile_size + self.metatile_size = tile_size * metatiling + # size in map units + self.x_size = float(round(self.right - self.left, ROUND)) + self.y_size = float(round(self.top - self.bottom, ROUND)) + + @property + def type(self): + warnings.warn(DeprecationWarning("'type' attribute is deprecated")) + return self.grid.type + + @property + def srid(self): + warnings.warn(DeprecationWarning("'srid' attribute is deprecated")) + return self.grid.srid + + def tile(self, zoom=None, row=None, col=None): + """ + Return Tile object of this TilePyramid. + + - zoom: zoom level + - row: tile matrix row + - col: tile matrix column + """ + return Tile(self, zoom, row, col) + + def matrix_width(self, zoom): + """ + Tile matrix width (number of columns) at zoom level. + + - zoom: zoom level + """ + validate_zoom(zoom) + width = int(math.ceil(self.grid.shape.width * 2 ** (zoom) / self.metatiling)) + return 1 if width < 1 else width + + def matrix_height(self, zoom): + """ + Tile matrix height (number of rows) at zoom level. + + - zoom: zoom level + """ + validate_zoom(zoom) + height = int(math.ceil(self.grid.shape.height * 2 ** (zoom) / self.metatiling)) + return 1 if height < 1 else height + + def tile_x_size(self, zoom): + """ + Width of a tile in SRID units at zoom level. + + - zoom: zoom level + """ + warnings.warn(DeprecationWarning("tile_x_size is deprecated")) + validate_zoom(zoom) + return round(self.x_size / self.matrix_width(zoom), ROUND) + + def tile_y_size(self, zoom): + """ + Height of a tile in SRID units at zoom level. + + - zoom: zoom level + """ + warnings.warn(DeprecationWarning("tile_y_size is deprecated")) + validate_zoom(zoom) + return round(self.y_size / self.matrix_height(zoom), ROUND) + + def tile_width(self, zoom): + """ + Tile width in pixel. + + - zoom: zoom level + """ + warnings.warn(DeprecationWarning("tile_width is deprecated")) + validate_zoom(zoom) + matrix_pixel = 2 ** (zoom) * self.tile_size * self.grid.shape.width + tile_pixel = self.tile_size * self.metatiling + return matrix_pixel if tile_pixel > matrix_pixel else tile_pixel + + def tile_height(self, zoom): + """ + Tile height in pixel. + + - zoom: zoom level + """ + warnings.warn(DeprecationWarning("tile_height is deprecated")) + validate_zoom(zoom) + matrix_pixel = 2 ** (zoom) * self.tile_size * self.grid.shape.height + tile_pixel = self.tile_size * self.metatiling + return matrix_pixel if tile_pixel > matrix_pixel else tile_pixel + + def pixel_x_size(self, zoom): + """ + Width of a pixel in SRID units at zoom level. + + - zoom: zoom level + """ + validate_zoom(zoom) + return round( + (self.grid.right - self.grid.left) + / (self.grid.shape.width * 2**zoom * self.tile_size), + ROUND, + ) + + def pixel_y_size(self, zoom): + """ + Height of a pixel in SRID units at zoom level. + + - zoom: zoom level + """ + validate_zoom(zoom) + return round( + (self.grid.top - self.grid.bottom) + / (self.grid.shape.height * 2**zoom * self.tile_size), + ROUND, + ) + + def intersecting(self, tile): + """ + Return all tiles intersecting with tile. + + This helps translating between TilePyramids with different metatiling + settings. + + - tile: a Tile object + """ + return _tile_intersecting_tilepyramid(tile, self) + + def tiles_from_bounds(self, bounds=None, zoom=None, batch_by=None): + """ + Return all tiles intersecting with bounds. + + Bounds values will be cleaned if they cross the antimeridian or are + outside of the Northern or Southern tile pyramid bounds. + + - bounds: tuple of (left, bottom, right, top) bounding values in tile + pyramid CRS + - zoom: zoom level + - batch_by: yield tiles in row or column batches if activated + """ + validate_zoom(zoom) + if not isinstance(bounds, tuple) or len(bounds) != 4: + raise ValueError( + "bounds must be a tuple of left, bottom, right, top values" + ) + if not isinstance(bounds, Bounds): + bounds = Bounds(*bounds) + if self.is_global: + yield from _global_tiles_from_bounds(self, bounds, zoom, batch_by=batch_by) + else: + yield from _tiles_from_cleaned_bounds(self, bounds, zoom, batch_by=batch_by) + + def tiles_from_bbox(self, geometry=None, zoom=None, batch_by=None): + """ + All metatiles intersecting with given bounding box. + + - geometry: shapely geometry + - zoom: zoom level + """ + validate_zoom(zoom) + yield from self.tiles_from_bounds(geometry.bounds, zoom, batch_by=batch_by) + + def tiles_from_geom(self, geometry=None, zoom=None, batch_by=None, exact=False): + """ + Return all tiles intersecting with input geometry. + + - geometry: shapely geometry + - zoom: zoom level + """ + validate_zoom(zoom) + if geometry.is_empty: + return + if not geometry.is_valid: + raise ValueError("no valid geometry: %s" % geometry.type) + if geometry.geom_type == "Point": + if batch_by: + yield ( + self.tile_from_xy(geometry.x, geometry.y, zoom) for _ in range(1) + ) + else: + yield self.tile_from_xy(geometry.x, geometry.y, zoom) + elif geometry.geom_type in ( + "MultiPoint", + "LineString", + "MultiLineString", + "Polygon", + "MultiPolygon", + "GeometryCollection", + ): + if exact: + geometry = clip_geometry_to_srs_bounds(geometry, self) + if batch_by: + for batch in self.tiles_from_bbox( + geometry, zoom, batch_by=batch_by + ): + yield ( + tile + for tile in batch + if geometry.intersection(tile.bbox()).area + ) + else: + for tile in self.tiles_from_bbox(geometry, zoom, batch_by=batch_by): + if geometry.intersection(tile.bbox()).area: + yield tile + else: + prepared_geometry = prep(clip_geometry_to_srs_bounds(geometry, self)) + if batch_by: + for batch in self.tiles_from_bbox( + geometry, zoom, batch_by=batch_by + ): + yield ( + tile + for tile in batch + if prepared_geometry.intersects(tile.bbox()) + ) + else: + for tile in self.tiles_from_bbox(geometry, zoom, batch_by=batch_by): + if prepared_geometry.intersects(tile.bbox()): + yield tile + + def tile_from_xy(self, x=None, y=None, zoom=None, on_edge_use="rb"): + """ + Return tile covering a point defined by x and y values. + + - x: x coordinate + - y: y coordinate + - zoom: zoom level + - on_edge_use: determine which Tile to pick if Point hits a grid edge + - rb: right bottom (default) + - rt: right top + - lt: left top + - lb: left bottom + """ + validate_zoom(zoom) + if x < self.left or x > self.right or y < self.bottom or y > self.top: + raise ValueError("x or y are outside of grid bounds") + if on_edge_use not in ["lb", "rb", "rt", "lt"]: + raise ValueError("on_edge_use must be one of lb, rb, rt or lt") + return _tile_from_xy(self, x, y, zoom, on_edge_use=on_edge_use) + + def to_dict(self): + """ + Return dictionary representation of pyramid parameters. + """ + return dict( + grid=self.grid.to_dict(), + metatiling=self.metatiling, + tile_size=self.tile_size, + ) + + def from_dict(config_dict): + """ + Initialize TilePyramid from configuration dictionary. + """ + return TilePyramid(**config_dict) + + def __eq__(self, other): + return ( + isinstance(other, self.__class__) + and self.grid == other.grid + and self.tile_size == other.tile_size + and self.metatiling == other.metatiling + ) + + def __ne__(self, other): + return not self.__eq__(other) + + def __repr__(self): + return "TilePyramid(%s, tile_size=%s, metatiling=%s)" % ( + self.grid, + self.tile_size, + self.metatiling, + ) + + def __hash__(self): + return hash(repr(self) + repr(self.grid)) diff --git a/tilematrix/_types.py b/tilematrix/_types.py new file mode 100644 index 0000000..e10335b --- /dev/null +++ b/tilematrix/_types.py @@ -0,0 +1,5 @@ +from collections import namedtuple + +Bounds = namedtuple("Bounds", "left bottom right top") +Shape = namedtuple("Shape", "height width") +TileIndex = namedtuple("TileIndex", "zoom row col") diff --git a/tilematrix/tmx/__init__.py b/tilematrix/tmx/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tilematrix/tmx/main.py b/tilematrix/tmx/main.py new file mode 100644 index 0000000..2cb41af --- /dev/null +++ b/tilematrix/tmx/main.py @@ -0,0 +1,199 @@ +import click +import geojson +from shapely.geometry import box + +import tilematrix +from tilematrix import TilePyramid + + +@click.version_option(version=tilematrix.__version__, message="%(version)s") +@click.group() +@click.option( + "--pixelbuffer", + "-p", + nargs=1, + type=click.INT, + default=0, + help="Tile bounding box buffer in pixels. (default: 0)", +) +@click.option( + "--tile_size", + "-s", + nargs=1, + type=click.INT, + default=256, + help="Tile size in pixels. (default: 256)", +) +@click.option( + "--metatiling", + "-m", + nargs=1, + type=click.INT, + default=1, + help="TilePyramid metatile size. (default: 1)", +) +@click.option( + "--grid", + "-g", + type=click.Choice(["geodetic", "mercator"]), + default="geodetic", + help="TilePyramid base grid. (default: geodetic)", +) +@click.option( + "--output_format", + "-f", + type=click.Choice(["Tile", "WKT", "GeoJSON"]), + default="Tile", + help="Print Tile ID or Tile bounding box as WKT or GeoJSON. (default: Tile)", +) +@click.pass_context +def tmx(ctx, **kwargs): + ctx.obj = dict(**kwargs) + + +@tmx.command(short_help="Tile bounds.") +@click.argument("TILE", nargs=3, type=click.INT, required=True) +@click.pass_context +def bounds(ctx, tile): + """Print Tile bounds.""" + click.echo( + "%s %s %s %s" + % TilePyramid( + ctx.obj["grid"], + tile_size=ctx.obj["tile_size"], + metatiling=ctx.obj["metatiling"], + ) + .tile(*tile) + .bounds(pixelbuffer=ctx.obj["pixelbuffer"]) + ) + + +@tmx.command(short_help="Tile bounding box.") +@click.argument("TILE", nargs=3, type=click.INT, required=True) +@click.pass_context +def bbox(ctx, tile): + """Print Tile bounding box as geometry.""" + geom = ( + TilePyramid( + ctx.obj["grid"], + tile_size=ctx.obj["tile_size"], + metatiling=ctx.obj["metatiling"], + ) + .tile(*tile) + .bbox(pixelbuffer=ctx.obj["pixelbuffer"]) + ) + if ctx.obj["output_format"] in ["WKT", "Tile"]: + click.echo(geom) + elif ctx.obj["output_format"] == "GeoJSON": + click.echo(geojson.dumps(geom)) + + +@tmx.command(short_help="Tile from point.") +@click.argument("ZOOM", nargs=1, type=click.INT, required=True) +@click.argument("POINT", nargs=2, type=click.FLOAT, required=True) +@click.pass_context +def tile(ctx, point, zoom): + """Print Tile containing POINT..""" + tile = TilePyramid( + ctx.obj["grid"], + tile_size=ctx.obj["tile_size"], + metatiling=ctx.obj["metatiling"], + ).tile_from_xy(*point, zoom=zoom) + if ctx.obj["output_format"] == "Tile": + click.echo("%s %s %s" % tile.id) + elif ctx.obj["output_format"] == "WKT": + click.echo(tile.bbox(pixelbuffer=ctx.obj["pixelbuffer"])) + elif ctx.obj["output_format"] == "GeoJSON": + click.echo( + geojson.dumps( + geojson.FeatureCollection( + [ + geojson.Feature( + geometry=tile.bbox(pixelbuffer=ctx.obj["pixelbuffer"]), + properties=dict(zoom=tile.zoom, row=tile.row, col=tile.col), + ) + ] + ) + ) + ) + + +@tmx.command(short_help="Tiles from bounds.") +@click.argument("ZOOM", nargs=1, type=click.INT, required=True) +@click.argument("BOUNDS", nargs=4, type=click.FLOAT, required=True) +@click.pass_context +def tiles(ctx, bounds, zoom): + """Print Tiles from bounds.""" + tiles = TilePyramid( + ctx.obj["grid"], + tile_size=ctx.obj["tile_size"], + metatiling=ctx.obj["metatiling"], + ).tiles_from_bounds(bounds, zoom=zoom) + if ctx.obj["output_format"] == "Tile": + for tile in tiles: + click.echo("%s %s %s" % tile.id) + elif ctx.obj["output_format"] == "WKT": + for tile in tiles: + click.echo(tile.bbox(pixelbuffer=ctx.obj["pixelbuffer"])) + elif ctx.obj["output_format"] == "GeoJSON": + click.echo("{\n" ' "type": "FeatureCollection",\n' ' "features": [') + # print tiles as they come and only add comma if there is a next tile + try: + tile = next(tiles) + while True: + gj = " %s" % geojson.Feature( + geometry=tile.bbox(pixelbuffer=ctx.obj["pixelbuffer"]), + properties=dict(zoom=tile.zoom, row=tile.row, col=tile.col), + ) + try: + tile = next(tiles) + click.echo(gj + ",") + except StopIteration: + click.echo(gj) + raise + except StopIteration: + pass + click.echo(" ]\n" "}") + + +@tmx.command(short_help="Snap bounds to tile grid.") +@click.argument("ZOOM", nargs=1, type=click.INT, required=True) +@click.argument("BOUNDS", nargs=4, type=click.FLOAT, required=True) +@click.pass_context +def snap_bounds(ctx, bounds, zoom): + """nap bounds to tile grid.""" + click.echo( + "%s %s %s %s" + % tilematrix.snap_bounds( + bounds=bounds, + tile_pyramid=TilePyramid( + ctx.obj["grid"], + tile_size=ctx.obj["tile_size"], + metatiling=ctx.obj["metatiling"], + ), + zoom=zoom, + pixelbuffer=ctx.obj["pixelbuffer"], + ) + ) + + +@tmx.command(short_help="Snap bbox to tile grid.") +@click.argument("ZOOM", nargs=1, type=click.INT, required=True) +@click.argument("BOUNDS", nargs=4, type=click.FLOAT, required=True) +@click.pass_context +def snap_bbox(ctx, bounds, zoom): + """Snap bbox to tile grid.""" + click.echo( + box( + *tilematrix.snap_bounds( + bounds=bounds, + tile_pyramid=TilePyramid( + ctx.obj["grid"], + tile_size=ctx.obj["tile_size"], + metatiling=ctx.obj["metatiling"], + ), + zoom=zoom, + pixelbuffer=ctx.obj["pixelbuffer"], + ) + ) + ) diff --git a/tilepyramidcut/_funcs.py b/tilepyramidcut/_funcs.py new file mode 100644 index 0000000..0cf4d67 --- /dev/null +++ b/tilepyramidcut/_funcs.py @@ -0,0 +1,161 @@ +import math +from pyproj import Transformer +import json +from shapely.geometry import shape +from shapely import from_geojson +from shapely.geometry.polygon import Polygon +from tilematrix import Tile, TilePyramid + +geos = """{ + "id": "4187b9f5-0aee-497e-b21d-a67451ecc2c6", + "pk": "", + "name": "unsaved trip", + "type": "trip", + "access": "public", + "bboxPolygon": { + "type": "Feature", + "bbox": [ + null, + null, + null, + null + ], + "properties": {}, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + null, + null + ], + [ + null, + null + ], + [ + null, + null + ], + [ + null, + null + ], + [ + null, + null + ] + ] + ] + } + }, + "events": { + "pre": [], + "during": [], + "post": [] + }, + "featureCollection": { + "type": "FeatureCollection", + "features": [ + { + "id": "581486ab552bc6495318301f65570213", + "type": "Feature", + "properties": { + "isRoute": true, + "name": "route 1", + "color": "#ff0000", + "isVisible": true, + "length": 387.73835, + "colorOuter": "#FFFFFF", + "speed": 6, + "departureTime": "2023-12-26T17:57:21.000Z", + "arrivalTime": "2023-12-29T10:34:44.008Z" + }, + "geometry": { + "coordinates": [ + [ + -122.6223, + 48.39416 + ], + [ + -124.79846, + 48.49026 + ], + [ + -125.00316, + 47.39116 + ], + [ + -124.3922, + 46.79512 + ], + [ + -122.18454, + 46.89421 + ], + [ + -121.87906, + 47.85601 + ], + [ + -122.58451, + 48.34604 + ] + ], + "type": "LineString" + } + } + ], + "bbox": [ + null, + null, + null, + null + ] + }, + "routeCount": 1, + "routesTotalLength": 387.73835, + "routesTotalDurationInDays": 2.69263, + "markerCount": 0, + "trackCount": 0, + "tracksTotalLength": 0 +}""" + +geo = json.loads(geos) +fc = geo["featureCollection"] +feature = fc["features"][0] +polygon: Polygon = shape(feature["geometry"]) + +# kk = from_geojson(fc) + + + +lon = 95.99 +lat = 15.15 + +t4326 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) +utm = t4326.transform(lon, lat) +print (utm) + +t3857 = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True) +lonlat = t3857.transform(utm[0], utm[1]) +print (lonlat) + +xx = list(map(lambda p:p[0], polygon.coords)) +yy = list(map(lambda p:p[1], polygon.coords)) + +mm = t4326.transform(xx, yy) +meters = list(zip(mm[0], mm[1])) +polym = Polygon(meters) + +tp = TilePyramid("mercator") +tiles = tp.tiles_from_geom(polym, 5) +for tile in tiles: + print (tile) + + + + + + + diff --git a/tilepyramidcut/test.py b/tilepyramidcut/test.py new file mode 100644 index 0000000..437217d --- /dev/null +++ b/tilepyramidcut/test.py @@ -0,0 +1,31 @@ +"""Tile geometries and tiles from geometries.""" + +from types import GeneratorType + +# import pytest +from shapely.geometry import Point, Polygon, shape + +from tilematrix import Tile, TilePyramid + + +def test_top_left_coord(): + """Top left coordinate.""" + # tp = TilePyramid({ + # "shape": (1, 1), # tile rows and columns at zoom level 0 + # "bounds": (-180.0, -85.051129, 180.0, 85.051129), # pyramid bounds in pyramid CRS + # "is_global": True, # if false, no antimeridian handling + # "srs": {"epsg": 4326}, # EPSG code for CRS + # }) + + tp = TilePyramid("mercator") + + tile = tp.tile(5, 3, 3) + #tile = tp.tile(1, 0, 0) + print(tile.bounds()) + print(tile.left, tile.bottom, tile.right, tile.top ) + + tp = TilePyramid("geodetic") + tile = tp.tile(5, 3, 3) + assert tile.bounds() == (-163.125, 67.5, -157.5, 73.125) + +test_top_left_coord() \ No newline at end of file