#!/bin/zsh

###############################################
# CONFIG
###############################################

PYTHON="/opt/local/Library/Frameworks/Python.framework/Versions/3.11/bin/python3"

TERRAIN_ROOT="terrain_rgb"
TILES_ROOT="tiles"
LOG_DIR="logs/tiles"
RESUME_DIR="resume/tiles"

Z_MIN=5
Z_MAX=14

mkdir -p "$TILES_ROOT" "$LOG_DIR" "$RESUME_DIR"


###############################################
# CHECK PYTHON
###############################################

if [ ! -x "$PYTHON" ]; then
    echo "[ERRORE] Python non trovato: $PYTHON"
    exit 1
fi


###############################################
# PYTHON SCRIPT INLINE (NUOVO E CORRETTO)
###############################################

PY_SCRIPT="$(mktemp)"
cat > "$PY_SCRIPT" << 'EOF'
import os, sys, math
from concurrent.futures import ProcessPoolExecutor
import rasterio
from rasterio.warp import transform_bounds
from rasterio.enums import Resampling
from rasterio.windows import Window
from PIL import Image


# ---------- Tile math ----------
def latlon_to_tile(lat, lon, z):
    lat_rad = math.radians(lat)
    n = 2 ** z
    xtile = int((lon + 180) / 360 * n)
    ytile = int((1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n)
    return xtile, ytile

def tile_bounds_wgs84(z, x, y):
    n = 2 ** z
    lon_min = x / n * 360 - 180
    lon_max = (x + 1) / n * 360 - 180
    lat_min = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * (y + 1) / n))))
    lat_max = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * y / n))))
    return lon_min, lat_min, lon_max, lat_max


# ---------- Clamp window ----------
def clamp_window(window, src):
    row_off = max(0, int(window.row_off))
    col_off = max(0, int(window.col_off))
    height = max(0, int(min(window.height, src.height - row_off)))
    width  = max(0, int(min(window.width,  src.width - col_off)))
    if width <= 0 or height <= 0:
        return None
    return Window(col_off, row_off, width, height)


# ---------- Tile job ----------
def process_tile(args):
    tif_path, tiles_root, z, x, y = args
    try:
        with rasterio.open(tif_path) as src:
            # CRS check
            if src.crs is None:
                return False

            # Tile bounds in WGS84
            lon_min, lat_min, lon_max, lat_max = tile_bounds_wgs84(z, x, y)

            # Transform tile bounds  raster CRS
            left, bottom, right, top = transform_bounds(
                "EPSG:4326", src.crs, lon_min, lat_min, lon_max, lat_max, densify_pts=0
            )

            # Convert to pixel coords
            r0, c0 = src.index(left,  top)
            r1, c1 = src.index(right, bottom)

            row_min, row_max = sorted((r0, r1))
            col_min, col_max = sorted((c0, c1))

            window = clamp_window(
                Window(col_min, row_min, col_max - col_min, row_max - row_min),
                src
            )
            if window is None:
                return False

            # Read terrain RGB
            data = src.read(
                indexes=(1,2,3),
                out_shape=(3, 256, 256),
                window=window,
                resampling=Resampling.nearest
            )

            if data.size == 0:
                return False

            # Save tile
            out_dir = os.path.join(tiles_root, str(z), str(x))
            os.makedirs(out_dir, exist_ok=True)
            out_path = os.path.join(out_dir, f"{y}.png")

            Image.fromarray(data.transpose(1,2,0).astype("uint8"), "RGB").save(out_path)

        return True

    except Exception:
        return False


# ---------- Main ----------
def main():
    tif_path    = sys.argv[1]
    tiles_root  = sys.argv[2]
    z_min       = int(sys.argv[3])
    z_max       = int(sys.argv[4])
    max_workers = int(sys.argv[5])

    with rasterio.open(tif_path) as src:
        if src.crs is None:
            print("[WARN] CRS mancante:", tif_path)
            return

        left, bottom, right, top = transform_bounds(src.crs, "EPSG:4326", *src.bounds)

    jobs = []
    for z in range(z_min, z_max + 1):
        x_min, y_max = latlon_to_tile(bottom, left,  z)
        x_max, y_min = latlon_to_tile(top,    right, z)

        x0, x1 = sorted((x_min, x_max))
        y0, y1 = sorted((y_min, y_max))

        for x in range(x0, x1 + 1):
            for y in range(y0, y1 + 1):
                jobs.append((tif_path, tiles_root, z, x, y))

    if not jobs:
        print("[WARN] Nessuna tile da generare per", tif_path)
        return

    with ProcessPoolExecutor(max_workers=max_workers) as ex:
        for _ in ex.map(process_tile, jobs):
            pass


if __name__ == "__main__":
    main()
EOF


###############################################
# LOOP SU TUTTI I TERRAINRGB (SOLO DEM)
###############################################

for DEM_DIR in "$TERRAIN_ROOT"/*; do
    [ -d "$DEM_DIR" ] || continue

    BASENAME=$(basename "$DEM_DIR")

    # Processa SOLO DEM (salta EDM e FLM)
    case "$BASENAME" in
        *DEM) ;;  
        *) echo "[SKIP] Non  DEM: $BASENAME"; continue ;;
    esac

    TIF="$DEM_DIR/${BASENAME}_terrain_rgb.tif"

    if [ ! -f "$TIF" ]; then
        echo "[WARN] Nessun TIFF TerrainRGB in $DEM_DIR"
        continue
    fi

    OUT_DIR="$TILES_ROOT/$BASENAME"
    mkdir -p "$OUT_DIR"

    RESUME_FILE="$RESUME_DIR/${BASENAME}.done"
    LOGFILE="$LOG_DIR/${BASENAME}.log"

    if [ -f "$RESUME_FILE" ]; then
        echo "[SKIP] Tiles gi generati per $BASENAME"
        continue
    fi

    echo "[INFO] Genero tiles per $BASENAME"

    {
        echo "----------------------------------------"
        echo "DEM: $BASENAME"
        echo "Input: $TIF"
        echo "Output: $OUT_DIR"
        echo "Zoom: $Z_MIN$Z_MAX"
        echo "Start: $(date)"
        echo "----------------------------------------"

        CORES=$(sysctl -n hw.ncpu 2>/dev/null || echo 4)

        "$PYTHON" "$PY_SCRIPT" "$TIF" "$OUT_DIR" "$Z_MIN" "$Z_MAX" "$CORES"
        STATUS=$?

        echo "End: $(date)"
        echo "Status: $STATUS"
        echo "----------------------------------------"

    } >> "$LOGFILE" 2>&1


    if [ "$STATUS" -ne 0 ]; then
        echo "[ERRORE] Tiles falliti per $BASENAME"
        continue
    fi

    touch "$RESUME_FILE"
    echo "[OK] Tiles creati per $BASENAME"

done

rm -f "$PY_SCRIPT"
echo "[INFO] Fase 2 completata"