Tiles prepare.md

Предобработка массива тайлов

Пояснения к данной статье.

Получение полигона области интереса

Используется формат GeoJSON. Требуемую область можно вручную выделить во множество онлайн редакторов или в одной из ГИС и экспортировать. Так же можно получить необходимый полигон из OSM, если он там уже уже отмечен. Воспользуемся последним способом на примере одного из городов. Для этого я использую http://overpass-turbo.eu.

overpass-turbo

Границы городов в OSM обычно обозначаются мультиполигонами. Откроем на карте требуемую область и составим следующий запрос для города. В данном примере это будет Пятигорск:

relation["type"="multipolygon"]["name"="Пятигорск"]({{bbox}});
out geom;

После чего необходимо проверить результат и экспортировать его в GeoJSON.

В случаях, когда мультиполигон не используется, запрос к overpass выглядит так. Ещё можно добавть тег place=city для фильтрации.

way["name"="Кисловодск"]({{bbox}});
out geom;

Экспорт GeoJSON выпоняется так же.

Генерация списков тайлов

Теперь необходимо получить тайлы, которые входят в полученный полигон. Для начала преобразуем все координаты вершин полигона в номера тайлов, к которым они относятся при помощи скрипта на python. Для начала откроем GeoJSON и получим все вершины:

with open('export.geojson') as f:
  vertices = json.load(f)['features'][0]['geometry']['coordinates']

Преобразуем координаты вершин полигона в соответствующие им номера тайлов и найдём прямоугольник, в который вписан полигон:

for vertex in flatten(vertices):
    x, y = deg2num(vertex[1], vertex[0])
    vertices_tiles.append([x, y])
    if x < bbox["xmin"] or bbox["xmin"] == 0:
        bbox["xmin"] = x
    if y < bbox["ymin"] or bbox["ymin"] == 0:
        bbox["ymin"] = y
    if x > bbox["xmax"] or bbox["xmax"] == 0:
        bbox["xmax"] = x
    if y > bbox["ymax"] or bbox["ymax"] == 0:
        bbox["ymax"] = y

Для последующего преобразования зададим смещение номерам тайлов:

for tile in vertices_tiles:
    shifted_vertices_tiles.append((tile[0] - bbox["xmin"] + 1, tile[1] - bbox["ymin"] + 1))

Создадим пустой холст для заливки полигона:

img = Image.new("L", (bbox["xmax"] - bbox["xmin"] + 3, bbox["ymax"] - bbox["ymin"] + 3), 1)
draw = ImageDraw.Draw(img)

Зальём полигон и преобразуем изображение в массив:

draw.polygon(shifted_vertices_tiles, fill=255)
mask = numpy.array(img)

Сохраним изображения для проверки правильности заливки:

img.save("polygon.bmp")

Полученный результат:

Преобразуем обратно в номера тайлов. При использовании небольшого масштаба тайлы, которые полигон чуть захватывает, могут не попасть в заливку. Зальём так же такие спорные места.

for h in range(len(mask)-1):
    for w in range(len(mask[h])-1):
        if mask[h][w] == 255:
            tiles_in_polygon.append([w + bbox["xmin"] - 1, h + bbox["ymin"] - 1])
        else:
            if (mask[h][w+1] == 255 and mask[h+1][w] == 255) or (mask[h][w-1] == 255 and mask[h+1][w] == 255) or (mask[h-1][w] == 255 and mask[h][w-1] == 255) or (mask[h][w+1] == 255 and mask[h-1][w] == 255):
                tiles_in_polygon.append([w + bbox["xmin"] - 1, h + bbox["ymin"] - 1])

Сохраним ссылки на тайлы:

with open("tiles.txt", "w") as text_file:
    for tile in tiles_in_polygon:
        text_file.write("{}/{}/{}/{}.png \r\n".format(tiles_server, zoom, tile[0], tile[1]))

Скрипт полностью.

Tiles prepare.py

import json
import math
import numpy
from collections import Iterable
from PIL import Image, ImageDraw, ImageFilter, ImagePath 

zoom = 20 # Требуемый масштаб
tiles_server = "https://tile.openstreetmap.org" # Адрес тайлового сервера

vertices_tiles = []
shifted_vertices_tiles = []
tiles_in_polygon = []
bbox = {
    "xmin": 0,
    "ymin": 0,
    "xmax": 0,
    "ymax": 0
}

# Взято с OSM Wiki
def deg2num(lat_deg, lon_deg):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** 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 flatten(lis):
    for item in lis:
        if isinstance(item, (list, dict, tuple)) and len(item) == 2:
            yield item
        else:
            yield from flatten(item)

with open('export.geojson') as f:
  vertices = json.load(f)['features'][0]['geometry']['coordinates']

for vertex in flatten(vertices):
    x, y = deg2num(vertex[1], vertex[0])
    vertices_tiles.append([x, y])
    if x < bbox["xmin"] or bbox["xmin"] == 0:
        bbox["xmin"] = x
    if y < bbox["ymin"] or bbox["ymin"] == 0:
        bbox["ymin"] = y
    if x > bbox["xmax"] or bbox["xmax"] == 0:
        bbox["xmax"] = x
    if y > bbox["ymax"] or bbox["ymax"] == 0:
        bbox["ymax"] = y

for tile in vertices_tiles:
    shifted_vertices_tiles.append((tile[0] - bbox["xmin"] + 1, tile[1] - bbox["ymin"] + 1))

img = Image.new("L", (bbox["xmax"] - bbox["xmin"] + 3, bbox["ymax"] - bbox["ymin"] + 3), 1)
draw = ImageDraw.Draw(img)

draw.polygon(shifted_vertices_tiles, fill=255)
mask = numpy.array(img)

img.save("polygon.bmp")

for h in range(len(mask)-1):
    for w in range(len(mask[h])-1):
        if mask[h][w] == 255:
            tiles_in_polygon.append([w + bbox["xmin"] - 1, h + bbox["ymin"] - 1])
        else:
            if (mask[h][w+1] == 255 and mask[h+1][w] == 255) or (mask[h][w-1] == 255 and mask[h+1][w] == 255) or (mask[h-1][w] == 255 and mask[h][w-1] == 255) or (mask[h][w+1] == 255 and mask[h-1][w] == 255):
                tiles_in_polygon.append([w + bbox["xmin"] - 1, h + bbox["ymin"] - 1])

with open("tiles.txt", "w") as text_file:
    for tile in tiles_in_polygon:
        text_file.write("{}/{}/{}/{}.png \r\n".format(tiles_server, zoom, tile[0], tile[1]))