smart-city-digital-twin-martinique/scripts/generate_martinique_mbtiles.py

#!/usr/bin/env python3
"""Generate mbtiles for Martinique from OpenStreetMap tiles.

Downloads raster tiles for Martinique area (zoom 0-14)
and packages them into an mbtiles file compatible with OpenRemote.

Usage: python3 generate_martinique_mbtiles.py

Requirements: pip install mbutil requests
Or: pip install sqlite3 (stdlib) + requests
"""

import sqlite3
import os
import sys
import time
import math
import requests
import gzip
from concurrent.futures import ThreadPoolExecutor, as_completed

# Martinique bounds
MIN_LON, MIN_LAT = -61.5, 14.0
MAX_LON, MAX_LAT = -60.5, 15.0

# Tile range
MIN_ZOOM = 0
MAX_ZOOM = 14

# Output file
OUTPUT = "/tmp/mapdata_martinique.mbtiles"

# Rate limiting for OSM tile usage policy
TILE_DELAY = 0.1  # 100ms between requests
OSM_TILE_URL = "https://tile.openstreetmap.org/{z}/{x}/{y}.png"
USER_AGENT = "SmartCityMartinique-MapGenerator/1.0"

def deg2num(lat_deg, lon_deg, zoom):
    """Convert lat/lon to tile numbers at given zoom level."""
    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)

def num2deg(xtile, ytile, zoom):
    """Convert tile numbers to lat/lon (top-left of tile)."""
    n = 1 << zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lat_deg, lon_deg)

def download_tile(z, x, y):
    """Download a single tile from OSM."""
    url = OSM_TILE_URL.format(z=z, x=x, y=y)
    try:
        resp = requests.get(url, headers={"User-Agent": USER_AGENT}, timeout=10)
        if resp.status_code == 200:
            return (z, x, y, resp.content)
        elif resp.status_code == 404:
            return None  # No tile here (ocean)
        else:
            print(f"  HTTP {resp.status_code} for {z}/{x}/{y}")
            return None
    except Exception as e:
        print(f"  Error {z}/{x}/{y}: {e}")
        return None

def create_mbtiles():
    """Create the mbtiles database."""
    if os.path.exists(OUTPUT):
        os.remove(OUTPUT)
    
    conn = sqlite3.connect(OUTPUT)
    c = conn.cursor()
    
    # Create tables
    c.execute("""CREATE TABLE map (
        zoom_level INTEGER,
        tile_column INTEGER,
        tile_row INTEGER,
        tile_id TEXT,
        grid_id TEXT
    )""")
    c.execute("""CREATE TABLE images (
        tile_data BLOB,
        tile_id TEXT
    )""")
    c.execute("""CREATE TABLE metadata (
        name TEXT,
        value TEXT
    )""")
    
    # Insert metadata
    metadata = [
        ("name", "Martinique OSM"),
        ("type", "overlay"),
        ("version", "1.0"),
        ("description", "OpenStreetMap raster tiles for Martinique"),
        ("format", "png"),
        ("tilejson", "2.1.0"),
        ("scheme", "xyz"),
        ("minzoom", str(MIN_ZOOM)),
        ("maxzoom", str(MAX_ZOOM)),
        ("bounds", f"{MIN_LON},{MIN_LAT},{MAX_LON},{MAX_LAT}"),
        ("center", f"{(MIN_LON+MAX_LON)/2},{(MIN_LAT+MAX_LAT)/2},10"),
        ("attribution", "&copy; <a href='https://www.openstreetmap.org/copyright'>OpenStreetMap</a> contributors"),
    ]
    c.executemany("INSERT INTO metadata VALUES (?, ?)", metadata)
    
    return conn, c

def main():
    print(f"Generating mbtiles for Martinique...")
    print(f"Bounds: {MIN_LON},{MIN_LAT} → {MAX_LON},{MAX_LAT}")
    print(f"Zoom: {MIN_ZOOM}-{MAX_ZOOM}")
    print(f"Output: {OUTPUT}")
    
    conn, c = create_mbtiles()
    
    # Calculate total tiles
    total_tiles = 0
    tiles_to_download = []
    for z in range(MIN_ZOOM, MAX_ZOOM + 1):
        x1, y1 = deg2num(MIN_LAT, MIN_LON, z)
        x2, y2 = deg2num(MAX_LAT, MAX_LON, z)
        x_min, x_max = min(x1, x2), max(x1, x2)
        y_min, y_max = min(y1, y2), max(y1, y2)
        
        # Add buffer of 1 tile
        x_min = max(0, x_min - 1)
        x_max = min((1 << z) - 1, x_max + 1)
        y_min = max(0, y_min - 1)
        y_max = min((1 << z) - 1, y_max + 1)
        
        for x in range(x_min, x_max + 1):
            for y in range(y_min, y_max + 1):
                tiles_to_download.append((z, x, y))
                total_tiles += 1
    
    print(f"Total tiles to download: {total_tiles}")
    
    # Download tiles (single-threaded for rate limiting)
    downloaded = 0
    failed = 0
    start = time.time()
    
    for i, (z, x, y) in enumerate(tiles_to_download):
        result = download_tile(z, x, y)
        if result:
            _, _, _, data = result
            tile_id = f"{z}_{x}_{y}"
            TMS_y = (1 << z) - 1 - y  # Convert XYZ to TMS
            
            c.execute("INSERT INTO images VALUES (?, ?)", (data, tile_id))
            c.execute("INSERT INTO map VALUES (?, ?, ?, ?, ?)", (z, x, TMS_y, tile_id, None))
            downloaded += 1
        else:
            # Store empty tile record
            tile_id = f"{z}_{x}_{y}"
            TMS_y = (1 << z) - 1 - y
            c.execute("INSERT INTO map VALUES (?, ?, ?, ?, ?)", (z, x, TMS_y, tile_id, None))
            failed += 1
        
        if (i + 1) % 50 == 0:
            conn.commit()
            elapsed = time.time() - start
            rate = (i + 1) / elapsed
            remaining = (total_tiles - i - 1) / rate if rate > 0 else 0
            print(f"  Progress: {i+1}/{total_tiles} ({downloaded} ok, {failed} skip) | {rate:.1f} tiles/s | ETA: {remaining:.0f}s")
        
        time.sleep(TILE_DELAY)
    
    conn.commit()
    elapsed = time.time() - start
    print(f"\nDone! {downloaded} tiles downloaded, {failed} skipped, in {elapsed:.0f}s")
    
    # Create index
    c.execute("CREATE INDEX idx_map ON map (zoom_level, tile_column, tile_row)")
    conn.commit()
    
    size_mb = os.path.getsize(OUTPUT) / (1024 * 1024)
    print(f"Output: {OUTPUT} ({size_mb:.1f} MB)")
    conn.close()

if __name__ == "__main__":
    main()