238 lines
6.2 KiB
Bash
Executable file
238 lines
6.2 KiB
Bash
Executable file
|
|
#!/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"
|
|
|