snippets/geo_tiles_assembler/plot_osm_map_with_coordinates.py

import matplotlib.pyplot as plt
import matplotlib.image as img
import requests
import numpy as np
import os
import math

def deg2num(lat_deg, lon_deg, zoom):
    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

#BOUNDS = ((49.17, 1.31), (49.18, 1.35))
BOUNDS = ((48.7, 2.2), (48.8, 2.3))
BOUNDS = ((48.8, 2.26), (48.82, 2.32))


DEFAULT_HEADERS = {"User-Agent": "Mozilla/5.0 (Windows NT 10.0; rv:127.0) Gecko/20100101 Firefox/127.0"}

iz = 0
for z in range(14, 15):
    # https://tile.openstreetmap.org/11/1029/700.png
    # 1029, 1030 going to the est, increasing the longitude angle (x)
    # 700, 701, going to the south, deacreasing the flatitude latitude angle (y)

    TILES_BOUNDS = (
        deg2num(BOUNDS[0][0], BOUNDS[0][1], z),
        deg2num(BOUNDS[1][0], BOUNDS[1][1], z)
    )
    print(TILES_BOUNDS)

    fig, ax = plt.subplots()

    maxx = 0
    maxy = 0
    xrange = range(TILES_BOUNDS[0][0], TILES_BOUNDS[1][0])
    yrange = range(TILES_BOUNDS[1][1], TILES_BOUNDS[0][1])
    print(xrange, yrange)

    iy = 0
    for y in reversed(yrange):
        ix = 0
        for x in xrange:
            # https://tile.openstreetmap.org/17/65956/44733.png
            url = f"https://tile.openstreetmap.org/{z}/{x}/{y}.png"
            dir_path = f'tiles/{z}/{x}'
            tile_path = dir_path + f'/{y}.png'
            img_data = 0
            print(tile_path)
            if not os.path.exists(tile_path):
                res = requests.get(url, headers=DEFAULT_HEADERS)
                print(res)
                if res.status_code != 200:
                    print(res.content)
                    exit()
                img_data = res.content
                os.makedirs(dir_path, exist_ok=True)
                with open(tile_path, 'wb') as handler:
                    handler.write(img_data)
            im = img.imread(tile_path)
            print("read")
            print(ix, iy)
            ax.imshow(im, aspect='auto', extent=(ix, ix+1, iy, iy+1), zorder=-1)
            print("imshow")
            ix += 1
        iy += 1
    ax.set_xlim(0, len(xrange))
    ax.set_ylim(0, len(yrange))
    ax.set_aspect('equal', adjustable='box')
    plt.show()
    break
    iz += 1