first commit

download_tiles.py

@@ -0,0 +1,163 @@
+#!/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()