Merge pull request #14 from eukarya-inc/feature/dev-11-tilelist

#11 implement tile_list modules

doc/plateauutils.rst

@@ -3,4 +3,5 @@ plateauutils パッケージ
 
 .. toctree::
 
- plateauutils.mesh_geocorder
+ plateauutils.mesh_geocorder
+ plateauutils.tile_list

doc/plateauutils.tile_list.rst

@@ -0,0 +1,18 @@
+plateauutils.tile_list パッケージ
+======================================
+
+plateauutils.tile_list.geo_to_tile モジュール
+--------------------------------------------------
+
+.. automodule:: plateauutils.tile_list.geo_to_tile
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+plateauutils.tile_list.polygon_to_tile_list モジュール
+---------------------------------------------------------------
+
+.. automodule:: plateauutils.tile_list.polygon_to_tile_list
+ :members:
+ :undoc-members:
+ :show-inheritance:

plateauutils/cli.py

@@ -12,6 +12,11 @@
  MeshException,
  MeshCodeException,
 )
+from plateauutils.tile_list.geo_to_tile import (
+ point_to_tile as _point_to_tile,
+ tile_to_polygon as _tile_to_polygon,
+ TileRangeException,
+)
 
 
 @click.group()
@@ -50,6 +55,35 @@ def meshcode_to_polygon(meshcode):
  click.echo("Error: {}".format(e))
 
 
+@cli.group()
+def tile_list():
+ """Tile list commands"""
+ pass
+
+
+@tile_list.command()
+@click.argument("longitude", required=True, type=float)
+@click.argument("latitude", required=True, type=float)
+@click.argument("zoom", required=True, type=int)
+@click.argument("ext", required=False, type=str, default=".mvt")
+def point_to_tile(longitude, latitude, zoom, ext):
+ """Convert a point to a tile"""
+ point = Point(longitude, latitude)
+ tile = _point_to_tile(point, zoom, ext)
+ click.echo(tile)
+
+
+@tile_list.command()
+@click.argument("tile_path", required=True, type=str)
+def tile_to_polygon(tile_path):
+ """Convert a tile to a polygon"""
+ try:
+ polygon = _tile_to_polygon(tile_path)
+ click.echo(polygon.wkt)
+ except TileRangeException as e:
+ click.echo("Error: {}".format(e))
+
+
 def main():
  cli()
 

plateauutils/tests/test_cli.py

@@ -1,5 +1,10 @@
 from click.testing import CliRunner
-from plateauutils.cli import point_to_meshcode, meshcode_to_polygon
+from plateauutils.cli import (
+ point_to_meshcode,
+ meshcode_to_polygon,
+ point_to_tile,
+ tile_to_polygon,
+)
 
 
 def test_point_to_meshcode():
@@ -31,3 +36,27 @@ def test_invalid_meshcode_to_polygon():
  result = runner.invoke(meshcode_to_polygon, ["53394547141"])
  assert result.exit_code == 0
  assert result.output == "Error: Mesh code must be 10 or less digits\n"
+
+
+def test_point_to_tile():
+ runner = CliRunner()
+ result = runner.invoke(point_to_tile, ["139.71475", "35.70078", "2"])
+ assert result.exit_code == 0
+ assert result.output == "2/3/1.mvt\n"
+
+
+def test_tile_to_polygon():
+ runner = CliRunner()
+ result = runner.invoke(tile_to_polygon, ["16/58211/25806.mvt"])
+ assert result.exit_code == 0
+ assert (
+ result.output
+ == "POLYGON ((139.7625732421875 35.68407153314097, 139.76806640625 35.68407153314097, 139.76806640625 35.679609609368576, 139.7625732421875 35.679609609368576, 139.7625732421875 35.68407153314097))\n"
+ )
+
+
+def test_invalid_tile_to_polygon():
+ runner = CliRunner()
+ result = runner.invoke(tile_to_polygon, ["1/2/1.mvt"])
+ assert result.exit_code == 0
+ assert result.output == "Error: Tile range must be less than 2\n"

plateauutils/tile_list/geo_to_tile.py

@@ -0,0 +1,91 @@
+import math
+from shapely.geometry import Point, Polygon
+
+
+def _deg2num(lat_deg, lon_deg, zoom):
+ # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
+ lat_rad = math.radians(lat_deg)
+ n = 1 << zoom
+ xtile = int((lon_deg + 180.0) / 360.0 * n)
+ ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
+ return xtile, ytile
+
+
+def _num2deg(xtile, ytile, zoom):
+ # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
+ n = 1 << zoom
+ lon_deg = xtile / n * 360.0 - 180.0
+ lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
+ lat_deg = math.degrees(lat_rad)
+ return lat_deg, lon_deg
+
+
+class TileRangeException(Exception):
+ pass
+
+
+def point_to_tile(point: Point, zoom: int, ext: str = ".mvt") -> str:
+ """緯度経度をタイルに変換して返す
+
+ Parameters
+ ----------
+ point : shapely.geometry.Point
+ 緯度経度を示すPointオブジェクト
+ zoom : int
+ ズームレベル
+ ext : str, optional
+ タイルの拡張子, by default ".mvt"
+
+ Returns
+ -------
+ str
+ タイルのパス
+ """
+ # 緯度経度の取得
+ longitude = point.x
+ latitude = point.y
+
+ # タイルの取得
+ xtile, ytile = _deg2num(latitude, longitude, zoom)
+
+ return f"{zoom}/{xtile}/{ytile}{ext}"
+
+
+def tile_to_polygon(tile_path: str) -> Polygon:
+ """
+ タイルのパスをポリゴンに変換して返す
+
+ Parameters
+ ----------
+ tile_path : str
+ タイルのパス
+
+ Returns
+ -------
+ shapely.geometry.Polygon
+ タイルのパスに対応するポリゴン
+ """
+
+ # タイルのパスからタイルの座標を取得
+ zoom, xtile, ytile = tile_path.split("/")[:3]
+ zoom = int(zoom)
+ xtile = int(xtile)
+ ytile = int(ytile.split(".")[0])
+ if zoom < 0 or xtile < 0 or ytile < 0:
+ raise TileRangeException("Tile range must be positive integer")
+ max_tile = 2**zoom
+ if xtile >= max_tile or ytile >= max_tile:
+ raise TileRangeException(f"Tile range must be less than {max_tile}")
+
+ # タイルの座標から緯度経度を取得
+ from_latitude, from_longitude = _num2deg(xtile, ytile, zoom)
+ to_latitude, to_longitude = _num2deg(xtile + 1, ytile + 1, zoom)
+
+ coords = (
+ (from_longitude, from_latitude),
+ (to_longitude, from_latitude),
+ (to_longitude, to_latitude),
+ (from_longitude, to_latitude),
+ (from_longitude, from_latitude),
+ )
+ return Polygon(coords)

plateauutils/tile_list/polygon_to_tile_list.py

@@ -0,0 +1,92 @@
+from itertools import product
+from plateauutils.abc.polygon_to_list import PolygonToList
+from plateauutils.tile_list.geo_to_tile import (
+ point_to_tile,
+ tile_to_polygon,
+)
+from shapely.geometry import Polygon
+
+
+class PolygonToTileList(PolygonToList):
+ """Polygonからタイルのリストを返すクラス
+
+ Parameters
+ ----------
+ polygon : shapely.geometry.Point
+ 対象となるポリゴン
+ zoom : int
+ ズームレベル
+ ext : str, optional
+ タイルの拡張子, by default ".mvt"
+ """
+
+ def __init__(self, polygon: Polygon, zoom: int, ext: str = ".mvt"):
+ super().__init__(polygon)
+ self.zoom = zoom
+ self.ext = ext
+ self.split()
+
+ def split(self):
+ """対象となるポリゴンを分割してタイルのリストを作成する"""
+ start_tile = point_to_tile(self.start_pos, self.zoom)
+ end_tile = point_to_tile(self.end_pos, self.zoom)
+ _, start_x, start_y = start_tile.split("/")
+ start_x = int(start_x)
+ start_y = int(start_y.split(".")[0])
+ _, end_x, end_y = end_tile.split("/")
+ end_x = int(end_x)
+ end_y = int(end_y.split(".")[0])
+ self.targets = _tile_list_range(
+ start_x, start_y, end_x, end_y, self.zoom, self.ext
+ )
+
+ def output(self) -> list:
+ """タイルのリストを出力する
+
+ Returns
+ -------
+ list
+ タイルのリスト
+ """
+ output_list = []
+ for i in self.targets:
+ polygon = tile_to_polygon(i)
+ # ポリゴンが交差しているかどうかを判定
+ if self.polygon.intersects(polygon):
+ output_list.append(i)
+ return sorted(output_list)
+
+
+def _tile_list_range(
+ start_x: int, start_y: int, end_x: int, end_y: int, zoom: int, ext: str
+) -> list:
+ """タイルのリストを作成する
+
+ Parameters
+ ----------
+ start_x : int
+ 開始x座標
+ start_y : int
+ 開始y座標
+ end_x : int
+ 終了x座標
+ end_y : int
+ 終了y座標
+ zoom : int
+ ズームレベル
+ ext : str
+ タイルの拡張子
+
+ Returns
+ -------
+ list
+ タイルのリスト
+ """
+ matrix = []
+ matrix.append(range(zoom, zoom + 1))
+ matrix.append(range(start_x, end_x + 1))
+ # y軸は上下逆(コンピュータ座標に合わせる)
+ matrix.append(range(end_y, start_y + 1))
+ combinations = product(*matrix)
+ tile_base = ["/".join(map(str, i)) for i in combinations]
+ return [i + ext for i in tile_base]

plateauutils/tile_list/tests/test_geo_to_tile.py

@@ -0,0 +1,47 @@
+import pytest
+from shapely.geometry import Polygon, Point
+
+
+def test_point_to_tile():
+ from plateauutils.tile_list.geo_to_tile import point_to_tile
+
+ point = Point(139.767125, 35.681236)
+ zoom = 16
+ tile_path = point_to_tile(point, zoom)
+ assert tile_path == "16/58211/25806.mvt"
+
+ zoom = 0
+ tile_path = point_to_tile(point, zoom)
+ assert tile_path == "0/0/0.mvt"
+
+
+def test_tile_to_polygon():
+ from plateauutils.tile_list.geo_to_tile import tile_to_polygon
+
+ tile_path = "16/58211/25806.mvt"
+ polygon = tile_to_polygon(tile_path)
+ assert (
+ polygon.wkt
+ == "POLYGON ((139.7625732421875 35.68407153314097, 139.76806640625 35.68407153314097, 139.76806640625 35.679609609368576, 139.7625732421875 35.679609609368576, 139.7625732421875 35.68407153314097))"
+ )
+
+ tile_path = "0/0/0.mvt"
+ polygon = tile_to_polygon(tile_path)
+ assert (
+ polygon.wkt
+ == "POLYGON ((-180 85.0511287798066, 180 85.0511287798066, 180 -85.0511287798066, -180 -85.0511287798066, -180 85.0511287798066))"
+ )
+
+
+def test_invalid_tile_to_polygon():
+ from plateauutils.tile_list.geo_to_tile import tile_to_polygon, TileRangeException
+
+ tile_path = "-1/0/0.mvt"
+ with pytest.raises(TileRangeException) as e:
+ tile_to_polygon(tile_path)
+ assert str(e.value) == "Tile range must be positive integer"
+
+ tile_path = "1/2/0.mvt"
+ with pytest.raises(TileRangeException) as e:
+ tile_to_polygon(tile_path)
+ assert str(e.value) == "Tile range must be less than 2"