#!/usr/bin/env python3
"""Fetch + stitch Esri World Imagery tiles (no API key) around a lat/lon.

Esri World Imagery over Japan is high-res aerial. Tile scheme is standard
Web Mercator XYZ, same as OSM.
"""
import math, sys, urllib.request, io
from PIL import Image

URL = "https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
UA = {"User-Agent": "Mozilla/5.0 (X11; Linux x86_64) research/1.0"}


def deg2tile(lat, lon, z):
    lat_r = math.radians(lat)
    n = 2.0 ** z
    xt = (lon + 180.0) / 360.0 * n
    yt = (1.0 - math.asinh(math.tan(lat_r)) / math.pi) / 2.0 * n
    return xt, yt


def fetch(z, x, y):
    u = URL.format(z=z, x=x, y=y)
    req = urllib.request.Request(u, headers=UA)
    with urllib.request.urlopen(req, timeout=30) as r:
        return Image.open(io.BytesIO(r.read())).convert("RGB")


def build(lat, lon, z, half, out):
    xt, yt = deg2tile(lat, lon, z)
    cx, cy = int(xt), int(yt)
    n = half * 2 + 1
    canvas = Image.new("RGB", (256 * n, 256 * n))
    ok = 0
    for i, x in enumerate(range(cx - half, cx + half + 1)):
        for j, y in enumerate(range(cy - half, cy + half + 1)):
            try:
                canvas.paste(fetch(z, x, y), (i * 256, j * 256))
                ok += 1
            except Exception as e:
                print(f"  tile {z}/{x}/{y} FAIL {e}", file=sys.stderr)
    canvas.save(out, quality=92)
    # metres per pixel at this latitude
    mpp = 156543.03392 * math.cos(math.radians(lat)) / (2 ** z)
    span = mpp * 256 * n
    print(f"{out}  z={z} tiles={ok}/{n*n} size={canvas.size} "
          f"mpp={mpp:.3f} span={span:.0f}m")
    return canvas


if __name__ == "__main__":
    LAT, LON = 35.46246657, 139.68149001
    build(LAT, LON, 17, 2, "/tmp/daikoku_z17.jpg")   # wide: whole loop
    build(LAT, LON, 18, 2, "/tmp/daikoku_z18.jpg")   # closer
    build(LAT, LON, 19, 2, "/tmp/daikoku_z19.jpg")   # detail: stalls