su2-img/scripts/download_tiles.py

#!/usr/bin/env python3
"""
Download satellite tiles from Mapy.com REST API for Hradec Králové.
Zoom 18 ≈ 0.38 m/pixel at lat 50° — good for car detection.

API key: https://developer.mapy.com/  (zdarma po registraci)
"""

import math
import os
import time
import argparse
from pathlib import Path
from concurrent.futures import ThreadPoolExecutor, as_completed

import requests
from PIL import Image

# Hradec Králové bounding box (min_lon, min_lat, max_lon, max_lat)
HK_BBOX = (15.7800, 50.1800, 15.8800, 50.2400)

# Mapy.com REST API v1 — mapset "aerial" = satelitní snímky
# Aerial: ČR až zoom 20, zbytek světa jen zoom 13
TILE_URL = "https://api.mapy.com/v1/maptiles/aerial/256/{z}/{x}/{y}?lang=cs&apikey={apikey}"
TILE_SIZE = 256
DEFAULT_ZOOM = 18
OUT_DIR = Path("tiles")


def lat_lon_to_tile(lat: float, lon: float, zoom: int) -> tuple[int, int]:
    """Convert lat/lon to tile x, y at given zoom (Slippy Map convention)."""
    n = 2 ** zoom
    x = int((lon + 180) / 360 * n)
    lat_rad = math.radians(lat)
    y = int((1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n)
    return x, y


def tile_to_lat_lon(x: int, y: int, zoom: int) -> tuple[float, float]:
    """Return NW corner (lat, lon) of tile."""
    n = 2 ** zoom
    lon = x / n * 360 - 180
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    lat = math.degrees(lat_rad)
    return lat, lon


def ground_resolution(lat: float, zoom: int) -> float:
    """Return meters per pixel at given lat and zoom."""
    return 156543.03 * math.cos(math.radians(lat)) / (2 ** zoom)


def download_tile(x: int, y: int, zoom: int, out_dir: Path, session: requests.Session, apikey: str) -> Path | None:
    path = out_dir / f"{zoom}_{x}_{y}.jpg"
    if path.exists():
        return path

    url = TILE_URL.format(z=zoom, x=x, y=y, apikey=apikey)
    try:
        r = session.get(url, timeout=15)
        r.raise_for_status()
        path.write_bytes(r.content)
        return path
    except Exception as e:
        print(f"  FAIL {x},{y}: {e}")
        return None


def stitch_tiles(tiles: dict[tuple, Path], x_range: range, y_range: range, zoom: int, out_path: Path):
    """Stitch downloaded tiles into one big image."""
    cols = len(x_range)
    rows = len(y_range)
    canvas = Image.new("RGB", (cols * TILE_SIZE, rows * TILE_SIZE))

    for j, y in enumerate(y_range):
        for i, x in enumerate(x_range):
            key = (x, y)
            if key in tiles and tiles[key]:
                try:
                    img = Image.open(tiles[key])
                    canvas.paste(img, (i * TILE_SIZE, j * TILE_SIZE))
                except Exception as e:
                    print(f"  Stitch error {x},{y}: {e}")

    canvas.save(out_path)
    print(f"Stitched image saved: {out_path}  ({canvas.width}x{canvas.height} px)")
    return canvas


def main():
    parser = argparse.ArgumentParser(description="Download Mapy.cz satellite tiles for Hradec Králové")
    parser.add_argument("--zoom", type=int, default=DEFAULT_ZOOM,
                        help=f"Zoom level (default: {DEFAULT_ZOOM})")
    parser.add_argument("--bbox", nargs=4, type=float,
                        metavar=("MIN_LON", "MIN_LAT", "MAX_LON", "MAX_LAT"),
                        default=HK_BBOX,
                        help="Bounding box (default: Hradec Králové)")
    parser.add_argument("--out", type=Path, default=OUT_DIR,
                        help="Output directory for tiles")
    parser.add_argument("--stitch", action="store_true",
                        help="Stitch all tiles into one image")
    parser.add_argument("--workers", type=int, default=8,
                        help="Parallel download workers")
    parser.add_argument("--delay", type=float, default=0.05,
                        help="Delay between requests (seconds)")
    parser.add_argument("--apikey", type=str, default=os.environ.get("MAPY_API_KEY", ""),
                        help="Mapy.com API key (nebo env MAPY_API_KEY)")
    args = parser.parse_args()

    if not args.apikey:
        parser.error("API klíč je povinný — použij --apikey nebo nastav env MAPY_API_KEY")

    min_lon, min_lat, max_lon, max_lat = args.bbox
    zoom = args.zoom

    res = ground_resolution((min_lat + max_lat) / 2, zoom)
    print(f"Zoom {zoom} → ~{res:.2f} m/pixel")

    # Tile range (NW corner has smaller y value)
    x_min, y_max = lat_lon_to_tile(min_lat, min_lon, zoom)
    x_max, y_min = lat_lon_to_tile(max_lat, max_lon, zoom)

    x_range = range(x_min, x_max + 1)
    y_range = range(y_min, y_max + 1)
    total = len(x_range) * len(y_range)

    print(f"Tiles: {len(x_range)} cols × {len(y_range)} rows = {total} total")
    print(f"Area: ~{len(x_range)*TILE_SIZE*res/1000:.1f} km × {len(y_range)*TILE_SIZE*res/1000:.1f} km")

    args.out.mkdir(parents=True, exist_ok=True)

    session = requests.Session()
    session.headers.update({
        "User-Agent": "Mozilla/5.0 (compatible; satellite-downloader/1.0)",
        "Referer": "https://mapy.cz/",
    })

    tiles: dict[tuple, Path | None] = {}
    done = 0

    jobs = [(x, y) for y in y_range for x in x_range]

    with ThreadPoolExecutor(max_workers=args.workers) as pool:
        futures = {pool.submit(download_tile, x, y, zoom, args.out, session, args.apikey): (x, y)
                   for x, y in jobs}
        for future in as_completed(futures):
            x, y = futures[future]
            tiles[(x, y)] = future.result()
            done += 1
            if done % 20 == 0 or done == total:
                print(f"  {done}/{total} tiles downloaded", end="\r")
            time.sleep(args.delay)

    ok = sum(1 for v in tiles.values() if v)
    print(f"\nDownloaded: {ok}/{total} tiles → {args.out}/")

    if args.stitch:
        stitch_out = args.out / f"hk_z{zoom}_stitched.jpg"
        stitch_tiles(tiles, x_range, y_range, zoom, stitch_out)


if __name__ == "__main__":
    main()