From 96037355197d6f7b0327bcbbc85cb7a44237c5ee Mon Sep 17 00:00:00 2001 From: Fabio Micheluz Date: Thu, 12 Feb 2026 17:21:04 +0100 Subject: [PATCH] first commit --- README.md | 188 ++++++++++++++++++++++++++++++++++ bounds.py | 96 +++++++++++++++++ e.py | 86 ++++++++++++++++ i.py | 74 +++++++++++++ make_tiles_rio_tiler.py | 96 +++++++++++++++++ make_tiles_rio_tiler.py.old | 58 +++++++++++ make_tiles_rio_tiler.py.old2 | 71 +++++++++++++ mbtiles_set_metadata.sh | 139 +++++++++++++++++++++++++ rgb2alt.py | 19 ++++ rgbify_batch.py | 49 +++++++++ set_mbtiles_metadata.py | 130 +++++++++++++++++++++++ set_mbtiles_metadata_dir.py | 118 +++++++++++++++++++++ set_mbtiles_metadata_multi.py | 85 +++++++++++++++ verify_and_fix_mbtiles.py | 144 ++++++++++++++++++++++++++ verify_bounds_tileservergl.py | 122 ++++++++++++++++++++++ verify_mbtiles_metadata.py | 122 ++++++++++++++++++++++ w.py | 68 ++++++++++++ 17 files changed, 1665 insertions(+) create mode 100644 README.md create mode 100644 bounds.py create mode 100644 e.py create mode 100644 i.py create mode 100644 make_tiles_rio_tiler.py create mode 100644 make_tiles_rio_tiler.py.old create mode 100644 make_tiles_rio_tiler.py.old2 create mode 100755 mbtiles_set_metadata.sh create mode 100644 rgb2alt.py create mode 100644 rgbify_batch.py create mode 100644 set_mbtiles_metadata.py create mode 100644 set_mbtiles_metadata_dir.py create mode 100644 set_mbtiles_metadata_multi.py create mode 100644 verify_and_fix_mbtiles.py create mode 100644 verify_bounds_tileservergl.py create mode 100644 verify_mbtiles_metadata.py create mode 100644 w.py diff --git a/README.md b/README.md new file mode 100644 index 0000000..16e192c --- /dev/null +++ b/README.md @@ -0,0 +1,188 @@ +# Creazione MbTiles per tileserver GL utilizzando Copernicus glo30 + +## Download dei file + +πŸš€ Scaricare i dati da mirror Copernicus GLO‑30 +``` +aws s3 sync s3://copernicus-dem-30m ./glo30 --no-sign-request +``` +Crea glo30 dir con tutto il mondo + +ma non servono tutti i dati quindi questi sono i tre programmi python per +scaricare italia europa e mondo in base a quello che vuoi ottenere + +``` +python3 i.py # italia +python3 e.py # europa +python3 w.py # mondo +``` + +## Trasformare i DEM geotiff in tiff RGB + +con la funzione rgbify_batch.py puoi creare le versioni tif in formato RGB + +``` +python3 rgbify_batch.py glo30_italia glo30_italia_rgb +``` + +il comando per la singola tiff sarebbe + +``` +rio rgbify -b -10000 -i 0.1 a.tif a_rgb.tif +``` + +## Trasformare gli RGB in tiles + +corretto +``` +python3 make_tiles_rio_tiler.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif tiles_x +``` + +## Trasformare i tiles in mbtiles + +``` +mb-util --scheme=xyz --image_format=png tiles_x/ x.mbtiles +``` + +## Inserire i metadata in mbtiles + +``` +python3 set_mbtiles_metadata.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif xxx.mbtiles "Terrain RGB Adamello" +``` + +## Verificare + +``` +python3 verify_mbtiles_metadata.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif xxx.mbtiles +python3 verify_and_fix_mbtiles.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif xxx.mbtiles +python3 verify_bounds_tileservergl.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif xxx.mbtiles +``` + +## multipli + +multili con file +``` +python3 make_tiles_rio_tiler.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif tiles_m +python3 make_tiles_rio_tiler.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E010_00_DEM_rgb.tif tiles_m + +mb-util --scheme=xyz --image_format=png tiles_m/ m.mbtiles + +python3 set_mbtiles_metadata_multi.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E010_00_DEM_rgb.tif m.mbtiles +``` + + +multipli con dir +``` +python3 make_tiles_rio_tiler_multi.py dir_tif tiles_mm +mb-util --scheme=xyz --image_format=png tiles_mm/ mm.mbtiles +python3 set_mbtiles_metadata_dir.py dir_tif mm.mbtiles +``` + +ultima idea +``` +python3 make_tiles_rio_tiler.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E010_00_DEM_rgb.tif tiles_x +python3 make_tiles_rio_tiler.py glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif tiles_y +rsync -av tiles_x/ tiles_m/ +rsync -av --ignore-existing tiles_y/ tiles_m/ +mb-util --scheme=xyz --image_format=png tiles_m/ m.mbtiles +python3 set_mbtiles_metadata_dir.py dir_tif m.mbtiles +scp m.mbtiles orangepi@192.168.1.4:/home/nvme/dockerdata/tileserver1/terrain.mbtiles +``` + +## Controllo + + +Monte Adamello + +Latitude: 46Β° 9' 19" N +Longitude: 10Β° 29' 47" E +Lat/Long (dec): 46.1554,10.4966 +KΓΆppen climate type: ET : Tundra +Elevation: 3,539m + +elevazione=-10000+(R*256^2+G*256+B)*0,1 + +estrarre l'altezza del monte dal geotiff +``` +gdallocationinfo -wgs84 -valonly glo30_italia/Copernicus_DSM_COG_10_N46_00_E010_00_DEM.tif 10.4966 46.1554 +``` +Γ¨ da come risultato + + 3496.90673828125 + +estrarre l'altezza del monte dal rgb +``` +gdallocationinfo -wgs84 glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E010_00_DEM_rgb.tif 10.4966 46.1554 +``` + +e da come report +``` +Report: + Location: (1788P,3041L) + Band 1: + Value: 2 + Band 2: + Value: 15 + Band 3: + Value: 57 +``` +utilizzando la funzione python che calcola altutudine da RGB si ottiene 3496.90 +``` +python3 rgb2alt.py 2 15 57 +``` + +verifichiamo con mbtiles + +``` +gdallocationinfo -wgs84 -valonly \ + -b 1 -b 2 -b 3 \ + aa.mbtiles 10.4966 46.1554 +``` + +``` +gdallocationinfo -wgs84 -valonly \ + -b 1 -b 2 -b 3 \ + x.mbtiles 10.4966 46.1554 +``` + +Punta Penia + +Dolomiti (3.343 m), sono approssimativamente: + +Latitudine: 46Β° 26' 03.84" N (46.43429) + +Longitudine: 11Β° 50' 58.58" E (11.85051) +``` +gdallocationinfo -wgs84 -valonly glo30_italia/Copernicus_DSM_COG_10_N46_00_E011_00_DEM.tif 11.85051 46.43429 +3236.89526367188 +gdallocationinfo -wgs84 glo30_italia_rgb/Copernicus_DSM_COG_10_N46_00_E011_00_DEM_rgb.tif 11.85051 46.43429 +Report: + Location: (3062P,2037L) + Band 1: + Value: 2 + Band 2: + Value: 5 + Band 3: + Value: 16 +python3 rgb2alt.py 2 5 16 +3236.800000000001 +gdallocationinfo -wgs84 -valonly -b 1 -b 2 -b 3 x.mbtiles 11.85051 46.43429 +2 +6 +68 +python3 rgb2alt.py 2 6 68 +3267.6 +``` + + +gdallocationinfo -wgs84 -valonly \ + -b 1 -b 2 -b 3 \ + aa.mbtiles 10.4966 46.1554 +2 +15 +147 +quindi 3505,9 + + + + diff --git a/bounds.py b/bounds.py new file mode 100644 index 0000000..f38cf7d --- /dev/null +++ b/bounds.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python3 +import argparse, sqlite3, sys + +def fetch_one(cur, sql, args=()): + row = cur.execute(sql, args).fetchone() + return row[0] if row else None + +def main(): + ap = argparse.ArgumentParser( + description="Normalizza bounds (W,S,E,N) e imposta center/minzoom/maxzoom in un MBTiles." + ) + ap.add_argument("mbtiles", help="Percorso al file .mbtiles") + ap.add_argument("--zoom", type=int, default=None, + help="Zoom del center da impostare (default: usa center esistente o media min/max dalle tile)") + ap.add_argument("--minzoom", type=int, default=None, + help="Forza minzoom (default: dai livelli in 'tiles')") + ap.add_argument("--maxzoom", type=int, default=None, + help="Forza maxzoom (default: dai livelli in 'tiles')") + args = ap.parse_args() + + conn = sqlite3.connect(args.mbtiles) + cur = conn.cursor() + + # --- 1) leggi bounds correnti + b_raw = fetch_one(cur, "SELECT value FROM metadata WHERE name='bounds';") + if not b_raw: + sys.exit("ERRORE: 'bounds' mancante nella tabella metadata.") + + try: + x1, y1, x2, y2 = [float(t) for t in b_raw.split(",")] + except Exception: + sys.exit(f"ERRORE: bounds non parseable: {b_raw}") + + # normalizza ordine W,S,E,N (min/max su lon e lat) + W, E = sorted((x1, x2)) + S, N = sorted((y1, y2)) + bounds = f"{W:.6f},{S:.6f},{E:.6f},{N:.6f}" + + # --- 2) min/max zoom: da tiles salvo override + minz = fetch_one(cur, "SELECT MIN(zoom_level) FROM tiles;") + maxz = fetch_one(cur, "SELECT MAX(zoom_level) FROM tiles;") + + # fallback se tiles non esiste o vuota + try: + minz = int(minz) if minz is not None else None + maxz = int(maxz) if maxz is not None else None + except Exception: + minz = maxz = None + + if args.minzoom is not None: minz = int(args.minzoom) + if args.maxzoom is not None: maxz = int(args.maxzoom) + + # --- 3) center: lon/lat dal bbox; zoom = argomento > zoom esistente > media min/max > 0 + center_existing = fetch_one(cur, "SELECT value FROM metadata WHERE name='center';") + center_z_existing = None + if center_existing: + parts = [p.strip() for p in center_existing.split(",")] + if len(parts) >= 3 and parts[2] != "": + center_z_existing = parts[2] + + if args.zoom is not None: + z = int(args.zoom) + elif center_z_existing is not None: + try: + z = int(center_z_existing) + except Exception: + z = None + else: + z = int((minz + maxz) // 2) if (minz is not None and maxz is not None) else 0 + + ctr_lon = (W + E) / 2.0 + ctr_lat = (S + N) / 2.0 + center = f"{ctr_lon:.6f},{ctr_lat:.6f},{z}" + + # --- 4) scrittura: crea metadata se manca, REPLACE dei valori + cur.execute("CREATE TABLE IF NOT EXISTS metadata (name TEXT PRIMARY KEY, value TEXT)") + cur.execute("REPLACE INTO metadata(name,value) VALUES(?,?)", ("bounds", bounds)) + cur.execute("REPLACE INTO metadata(name,value) VALUES(?,?)", ("center", center)) + if minz is not None: + cur.execute("REPLACE INTO metadata(name,value) VALUES(?,?)", ("minzoom", str(minz))) + if maxz is not None: + cur.execute("REPLACE INTO metadata(name,value) VALUES(?,?)", ("maxzoom", str(maxz))) + conn.commit() + + # --- 5) stampa riepilogo + print("OK") + for k in ("bounds","center","minzoom","maxzoom"): + v = fetch_one(cur, "SELECT value FROM metadata WHERE name=?", (k,)) + print(f"{k} = {v}") + zminmax = cur.execute("SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;").fetchone() + if zminmax and all(v is not None for v in zminmax): + print(f"zoom (tiles) = {zminmax[0]}..{zminmax[1]}") + conn.close() + +if __name__ == "__main__": + main() diff --git a/e.py b/e.py new file mode 100644 index 0000000..7e9ddcd --- /dev/null +++ b/e.py @@ -0,0 +1,86 @@ +import boto3 +import os +import logging +from botocore import UNSIGNED +from botocore.config import Config + +BUCKET = "copernicus-dem-30m" +OUTPUT_DIR = "./glo30_europe" +LOG_FILE = "download_europe.log" + +# Europa approx: +lat_range = range(34, 73) # N34–N72 +lon_east = range(0, 46) # E000–E045 +lon_west = range(1, 26) # W001–W025 + +logging.basicConfig( + filename=LOG_FILE, + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", +) + +console = logging.getLogger("console") +console.setLevel(logging.INFO) +console.addHandler(logging.StreamHandler()) + +def log(msg): + logging.info(msg) + console.info(msg) + +os.makedirs(OUTPUT_DIR, exist_ok=True) + +s3 = boto3.client("s3", config=Config(signature_version=UNSIGNED)) + +def folder_name(lat, lon, east=True): + if east: + return f"Copernicus_DSM_COG_10_N{lat:02d}_00_E{lon:03d}_00_DEM/" + else: + return f"Copernicus_DSM_COG_10_N{lat:02d}_00_W{lon:03d}_00_DEM/" + +log("Inizio scansione directory europee...") + +found_folders = [] + +paginator = s3.get_paginator("list_objects_v2") + +for page in paginator.paginate(Bucket=BUCKET, Delimiter="/"): + for prefix in page.get("CommonPrefixes", []): + folder = prefix["Prefix"] + + # EASTERN EUROPE + for lat in lat_range: + for lon in lon_east: + if folder == folder_name(lat, lon, east=True): + found_folders.append(folder) + log(f"[FOUND] {folder}") + + # WESTERN EUROPE (longitudes W) + for lat in lat_range: + for lon in lon_west: + if folder == folder_name(lat, lon, east=False): + found_folders.append(folder) + log(f"[FOUND] {folder}") + +log(f"Trovate {len(found_folders)} cartelle europee") + +# ----------------------------- +# DOWNLOAD DEM PRINCIPALE +# ----------------------------- +for folder in found_folders: + tif_name = folder[:-1] + ".tif" + key = folder + tif_name.split("/")[-1] + out_path = os.path.join(OUTPUT_DIR, tif_name.split("/")[-1]) + + if os.path.exists(out_path): + log(f"[SKIP] {out_path} giΓ  presente") + continue + + log(f"[DOWNLOAD] {key}") + try: + s3.download_file(BUCKET, key, out_path) + log(f"[OK] {key} scaricato") + except Exception as e: + log(f"[ERROR] {key} β†’ {e}") + +log("Completato.") + diff --git a/i.py b/i.py new file mode 100644 index 0000000..5ea57d2 --- /dev/null +++ b/i.py @@ -0,0 +1,74 @@ +import boto3 +import os +import logging +from botocore import UNSIGNED +from botocore.config import Config + +BUCKET = "copernicus-dem-30m" +OUTPUT_DIR = "./glo30_italia" +LOG_FILE = "download.log" + +lat_range = range(36, 48) # N36–N47 +lon_range = range(6, 19) # E006–E018 + +logging.basicConfig( + filename=LOG_FILE, + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", +) + +console = logging.getLogger("console") +console.setLevel(logging.INFO) +console.addHandler(logging.StreamHandler()) + +def log(msg): + logging.info(msg) + console.info(msg) + +os.makedirs(OUTPUT_DIR, exist_ok=True) + +s3 = boto3.client("s3", config=Config(signature_version=UNSIGNED)) + +def folder_name(lat, lon): + return f"Copernicus_DSM_COG_10_N{lat:02d}_00_E{lon:03d}_00_DEM/" + +log("Inizio scansione delle directory italiane...") + +found_folders = [] + +paginator = s3.get_paginator("list_objects_v2") + +for page in paginator.paginate(Bucket=BUCKET, Delimiter="/"): + for prefix in page.get("CommonPrefixes", []): + folder = prefix["Prefix"] + + # Controlla se la cartella Γ¨ italiana + for lat in lat_range: + for lon in lon_range: + if folder == folder_name(lat, lon): + found_folders.append(folder) + log(f"[FOUND] {folder}") + +log(f"Trovate {len(found_folders)} cartelle italiane") + +# ----------------------------- +# DOWNLOAD DEM PRINCIPALE +# ----------------------------- +for folder in found_folders: + tif_name = folder[:-1] + ".tif" # rimuove "/" e aggiunge .tif + key = folder + tif_name.split("/")[-1] + out_path = os.path.join(OUTPUT_DIR, tif_name.split("/")[-1]) + + if os.path.exists(out_path): + log(f"[SKIP] {out_path} giΓ  presente") + continue + + log(f"[DOWNLOAD] {key}") + try: + s3.download_file(BUCKET, key, out_path) + log(f"[OK] {key} scaricato") + except Exception as e: + log(f"[ERROR] {key} β†’ {e}") + +log("Completato.") + diff --git a/make_tiles_rio_tiler.py b/make_tiles_rio_tiler.py new file mode 100644 index 0000000..67a3e9d --- /dev/null +++ b/make_tiles_rio_tiler.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python3 +import os +import sys +import mercantile +import numpy as np +from rio_tiler.io import COGReader +from rio_tiler.utils import render + +# --------------------------------------------------------- +# USO: +# python3 make_tiles_rio_tiler.py input_rgb.tif tiles_out/ +# +# Versione finale: +# - BUFFER per evitare buchi tra quadranti +# - PADDING per generare tile anche se il TIFF copre solo parte della tile +# - Logging avanzato +# - Controllo tile vuote +# --------------------------------------------------------- + +ZOOM_MIN = 5 +ZOOM_MAX = 14 +BUFFER = 0.005 # espansione bounds per evitare tile mancanti +TILESIZE = 256 # tile standard +PADDING = 1 # genera tile anche se parziali + +def generate_tiles(input_tif, output_dir): + os.makedirs(output_dir, exist_ok=True) + + print(f"\n=== Generazione tile per: {input_tif} ===") + print(f"Zoom: {ZOOM_MIN}..{ZOOM_MAX}") + print(f"Buffer: {BUFFER}Β°") + print(f"Padding: {PADDING} pixel\n") + + with COGReader(input_tif) as cog: + + for z in range(ZOOM_MIN, ZOOM_MAX + 1): + print(f"\n== Zoom {z} ==") + + # Bounds originali del TIFF + bounds = cog.bounds + + # Espansione per evitare buchi tra quadranti + minx = bounds[0] - BUFFER + miny = bounds[1] - BUFFER + maxx = bounds[2] + BUFFER + maxy = bounds[3] + BUFFER + + # Tile da generare + tiles = list(mercantile.tiles(minx, miny, maxx, maxy, z)) + print(f"Tile da generare: {len(tiles)}") + + for t in tiles: + out_dir_zxy = os.path.join(output_dir, str(z), str(t.x)) + os.makedirs(out_dir_zxy, exist_ok=True) + out_path = os.path.join(out_dir_zxy, f"{t.y}.png") + + try: + tile, mask = cog.tile( + t.x, t.y, z, + tilesize=TILESIZE, + padding=PADDING, + resampling_method="nearest" + ) + except Exception as e: + print(f" ⚠️ Errore tile {z}/{t.x}/{t.y}: {e}") + continue + + # Controllo tile vuote (tutto zero o tutto trasparente) + if mask is not None and np.all(mask == 0): + print(f" ⚠️ Tile vuota saltata: {z}/{t.x}/{t.y}") + continue + + img = render(tile, mask=mask, img_format="PNG") + + with open(out_path, "wb") as f: + f.write(img) + + print(f" βœ“ {z}/{t.x}/{t.y}") + + print("\n=== Completato senza errori ===\n") + + +if __name__ == "__main__": + if len(sys.argv) != 3: + print("Uso: python3 make_tiles_rio_tiler.py ") + sys.exit(1) + + INPUT_TIF = sys.argv[1] + OUTPUT_DIR = sys.argv[2] + + if not os.path.isfile(INPUT_TIF): + print(f"Errore: file non trovato: {INPUT_TIF}") + sys.exit(1) + + generate_tiles(INPUT_TIF, OUTPUT_DIR) + diff --git a/make_tiles_rio_tiler.py.old b/make_tiles_rio_tiler.py.old new file mode 100644 index 0000000..ded6af6 --- /dev/null +++ b/make_tiles_rio_tiler.py.old @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 +import os +import sys +import mercantile +from rio_tiler.io import COGReader +from rio_tiler.utils import render + +# --------------------------------------------------------- +# USO: +# python3 make_tiles_rio_tiler.py input_rgb.tif tiles_out/ +# +# Genera tile PNG Terrain-RGB senza alterare i pixel. +# --------------------------------------------------------- + +ZOOM_MIN = 5 +ZOOM_MAX = 14 + +def generate_tiles(input_tif, output_dir): + os.makedirs(output_dir, exist_ok=True) + + with COGReader(input_tif) as cog: + for z in range(ZOOM_MIN, ZOOM_MAX + 1): + print(f"== Zoom {z} ==") + + bounds = cog.bounds + tiles = list(mercantile.tiles(bounds[0], bounds[1], bounds[2], bounds[3], z)) + + for t in tiles: + out_dir_zxy = os.path.join(output_dir, str(z), str(t.x)) + os.makedirs(out_dir_zxy, exist_ok=True) + out_path = os.path.join(out_dir_zxy, f"{t.y}.png") + + tile, mask = cog.tile(t.x, t.y, z, resampling_method="nearest") + + img = render(tile, mask=mask, img_format="PNG") + + with open(out_path, "wb") as f: + f.write(img) + + print(f"Tile {z}/{t.x}/{t.y} generata") + + print("Completato.") + + +if __name__ == "__main__": + if len(sys.argv) != 3: + print("Uso: python3 make_tiles_rio_tiler.py ") + sys.exit(1) + + INPUT_TIF = sys.argv[1] + OUTPUT_DIR = sys.argv[2] + + if not os.path.isfile(INPUT_TIF): + print(f"Errore: file non trovato: {INPUT_TIF}") + sys.exit(1) + + generate_tiles(INPUT_TIF, OUTPUT_DIR) + diff --git a/make_tiles_rio_tiler.py.old2 b/make_tiles_rio_tiler.py.old2 new file mode 100644 index 0000000..aaa5f00 --- /dev/null +++ b/make_tiles_rio_tiler.py.old2 @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 +import os +import sys +import mercantile +from rio_tiler.io import COGReader +from rio_tiler.utils import render + +# --------------------------------------------------------- +# USO: +# python3 make_tiles_rio_tiler.py input_rgb.tif tiles_out/ +# +# Genera tile PNG Terrain-RGB senza alterare i pixel. +# Corretto per evitare buchi tra quadranti Copernicus. +# --------------------------------------------------------- + +ZOOM_MIN = 5 +ZOOM_MAX = 14 +BUFFER = 0.001 # espansione per evitare tile mancanti ai bordi + +def generate_tiles(input_tif, output_dir): + os.makedirs(output_dir, exist_ok=True) + + with COGReader(input_tif) as cog: + for z in range(ZOOM_MIN, ZOOM_MAX + 1): + print(f"== Zoom {z} ==") + + # Bounds originali del TIFF + bounds = cog.bounds + + # Espansione per evitare buchi tra quadranti + minx = bounds[0] - BUFFER + miny = bounds[1] - BUFFER + maxx = bounds[2] + BUFFER + maxy = bounds[3] + BUFFER + + # Tile da generare + tiles = list(mercantile.tiles(minx, miny, maxx, maxy, z)) + + for t in tiles: + out_dir_zxy = os.path.join(output_dir, str(z), str(t.x)) + os.makedirs(out_dir_zxy, exist_ok=True) + out_path = os.path.join(out_dir_zxy, f"{t.y}.png") + + # Estrazione tile + tile, mask = cog.tile(t.x, t.y, z, resampling_method="nearest") + + # Render PNG + img = render(tile, mask=mask, img_format="PNG") + + with open(out_path, "wb") as f: + f.write(img) + + print(f"Tile {z}/{t.x}/{t.y} generata") + + print("Completato.") + + +if __name__ == "__main__": + if len(sys.argv) != 3: + print("Uso: python3 make_tiles_rio_tiler.py ") + sys.exit(1) + + INPUT_TIF = sys.argv[1] + OUTPUT_DIR = sys.argv[2] + + if not os.path.isfile(INPUT_TIF): + print(f"Errore: file non trovato: {INPUT_TIF}") + sys.exit(1) + + generate_tiles(INPUT_TIF, OUTPUT_DIR) + diff --git a/mbtiles_set_metadata.sh b/mbtiles_set_metadata.sh new file mode 100755 index 0000000..6442e2f --- /dev/null +++ b/mbtiles_set_metadata.sh @@ -0,0 +1,139 @@ +#!/usr/bin/env bash +set -euo pipefail + +# -------------------------------------------------------------------- +# Imposta metadati essenziali in un file .mbtiles (schema XYZ) +# - Calcola bounds (ovest,sud,est,nord) dal raster sorgente (GeoTIFF/COG) +# - Calcola center (lon,lat,zoom) +# - Legge min/max zoom dalla tabella 'tiles' e li scrive nella 'metadata' +# - Imposta name/format/type/attribution +# +# Requisiti: gdalinfo (GDAL), jq, sqlite3 +# -------------------------------------------------------------------- +# USO: +# ./mbtiles_set_metadata.sh \ +# --mbtiles aa.mbtiles \ +# --raster aa_rgb.tif \ +# [--name "Terrain RGB Adamello"] \ +# [--type overlay|baselayer] \ +# [--format png|jpg] \ +# [--zoom ] \ +# [--attribution " Copernicus DEM (GLO-30), ESA"] +# +# NOTE: +# - 'bounds' e 'center' sono calcolati in EPSG:4326 (lon/lat). +# - Se il raster non riassume esattamente lestensione delle tile, +# passa un raster "clip" coerente con larea tilata. +# - Lo script usa REPLACE INTO per compatibilit con SQLite < 3.24. +# -------------------------------------------------------------------- + +# ---------- parsing argomenti ---------- +MBTILES="" +RASTER="" +LAYER_NAME="" +LAYER_TYPE="overlay" +FORMAT="png" +CENTER_Z="" +ATTRIBUTION="" + +while [[ $# -gt 0 ]]; do + case "$1" in + --mbtiles) MBTILES="$2"; shift 2;; + --raster) RASTER="$2"; shift 2;; + --name) LAYER_NAME="$2"; shift 2;; + --type) LAYER_TYPE="$2"; shift 2;; + --format) FORMAT="$2"; shift 2;; + --zoom) CENTER_Z="$2"; shift 2;; + --attribution) ATTRIBUTION="$2"; shift 2;; + -h|--help) + sed -n '1,80p' "$0"; exit 0;; + *) + echo "Argomento sconosciuto: $1" >&2; exit 1;; + esac +done + +if [[ -z "$MBTILES" || -z "$RASTER" ]]; then + echo "Uso: $0 --mbtiles --raster [opzioni]" >&2 + exit 1 +fi + +if [[ ! -f "$MBTILES" ]]; then + echo "File MBTiles non trovato: $MBTILES" >&2; exit 1 +fi +if [[ ! -f "$RASTER" ]]; then + echo "Raster non trovato: $RASTER" >&2; exit 1 +fi + +command -v gdalinfo >/dev/null || { echo "gdalinfo non trovato"; exit 1; } +command -v jq >/dev/null || { echo "jq non trovato"; exit 1; } +command -v sqlite3 >/dev/null || { echo "sqlite3 non trovato"; exit 1; } + +# ---------- step 1: calcolo bounds dal raster (EPSG:4326) ---------- +# Usiamo wgs84Extent se presente, altrimenti cornerCoordinates. +# gdalinfo -json fornisce geometrie e corner coords in lon/lat. +export LC_ALL=C +export LANG=C + +GDAL_JSON="$(gdalinfo -json "$RASTER")" + +has_wgs84=$(echo "$GDAL_JSON" | jq 'has("wgs84Extent")') +if [[ "$has_wgs84" == "true" ]]; then + # Polygon ring: [ [minLon,minLat], [maxLon,minLat], [maxLon,maxLat], [minLon,maxLat], ... ] + read -r MINLON MINLAT MAXLON MAXLAT <<<"$(echo "$GDAL_JSON" \ + | jq -r '.wgs84Extent.coordinates[0] | + [.[0][0], .[0][1], .[2][0], .[2][1]] | @tsv')" +else + # cornerCoordinates: upperLeft, lowerLeft, upperRight, lowerRight (lon,lat) + read -r MINLON MINLAT MAXLON MAXLAT <<<"$(echo "$GDAL_JSON" \ + | jq -r '.cornerCoordinates as $c | + [$c.lowerLeft[0], $c.lowerLeft[1], $c.upperRight[0], $c.upperRight[1]] | @tsv')" +fi + +# Arrotonda a 6 decimali per compattare la stringa +printf -v BOUNDS "%.6f,%.6f,%.6f,%.6f" "$MINLON" "$MINLAT" "$MAXLON" "$MAXLAT" + +# ---------- step 2: calcolo center (lon,lat,zoom) ---------- +# lon,lat = met dei bounds; zoom: se non fornito, media tra min e max zoom presenti nelle tile +CTR_LON=$(awk -v a="$MINLON" -v b="$MAXLON" 'BEGIN{printf "%.6f",(a+b)/2.0}') +CTR_LAT=$(awk -v a="$MINLAT" -v b="$MAXLAT" 'BEGIN{printf "%.6f",(a+b)/2.0}') + +# Leggi min/max zoom dalla tabella 'tiles' +read -r TMIN TMAX <<<"$(sqlite3 "$MBTILES" "SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;")" +if [[ -z "$TMIN" || -z "$TMAX" ]]; then + echo "Attenzione: tabella 'tiles' vuota o mancante in $MBTILES" >&2 + TMIN=0; TMAX=0 +fi + +if [[ -z "${CENTER_Z}" ]]; then + CENTER_Z=$(( (TMIN + TMAX) / 2 )) +fi + +CENTER="${CTR_LON},${CTR_LAT},${CENTER_Z}" + +# ---------- step 3: definisci name/format/type/attribution ---------- +if [[ -z "$LAYER_NAME" ]]; then + # Usa il nome file come default + base="$(basename "$MBTILES")" + LAYER_NAME="${base%.*}" +fi + +# ---------- step 4: crea tabella metadata (se non esiste) e scrivi i metadati ---------- +sqlite3 "$MBTILES" < METADATA scritti in: $MBTILES" +sqlite3 "$MBTILES" "SELECT name,value FROM metadata ORDER BY name;" +echo +echo "Zoom effettivi nelle tile: $(sqlite3 "$MBTILES" "SELECT MIN(zoom_level)||'..'||MAX(zoom_level) FROM tiles;")" +echo "Bounds: $BOUNDS" +echo "Center: $CENTER" diff --git a/rgb2alt.py b/rgb2alt.py new file mode 100644 index 0000000..19edf7f --- /dev/null +++ b/rgb2alt.py @@ -0,0 +1,19 @@ +#!/usr/bin/env python3 + +def rgb_to_altitude(R, G, B, scale=0.1, offset=-10000): + value = (R * 65536) + (G * 256) + B + return value * scale + offset + +if __name__ == "__main__": + import sys + if len(sys.argv) != 4: + print("Uso: python3 rgb2alt.py R G B") + sys.exit(1) + + R = int(sys.argv[1]) + G = int(sys.argv[2]) + B = int(sys.argv[3]) + + alt = rgb_to_altitude(R, G, B) + print(alt) + diff --git a/rgbify_batch.py b/rgbify_batch.py new file mode 100644 index 0000000..8a9258e --- /dev/null +++ b/rgbify_batch.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python3 +import os +import sys +import subprocess + +def main(): + if len(sys.argv) != 3: + print("Uso: python3 rgbify_batch.py ") + sys.exit(1) + + input_dir = sys.argv[1] + output_dir = sys.argv[2] + + if not os.path.isdir(input_dir): + print(f"Errore: directory di input non trovata: {input_dir}") + sys.exit(1) + + os.makedirs(output_dir, exist_ok=True) + + # Scansiona tutti i .tif nella directory di input + for filename in os.listdir(input_dir): + if not filename.lower().endswith(".tif"): + continue + + in_path = os.path.join(input_dir, filename) + out_name = filename.replace(".tif", "_rgb.tif") + out_path = os.path.join(output_dir, out_name) + + print(f"[PROCESS] {filename} β†’ {out_name}") + + cmd = [ + "rio", "rgbify", + "-b", "-10000", + "-i", "0.1", + in_path, + out_path + ] + + try: + subprocess.run(cmd, check=True) + print(f"[OK] {out_name}") + except subprocess.CalledProcessError as e: + print(f"[ERROR] {filename}: {e}") + + print("Completato.") + +if __name__ == "__main__": + main() + diff --git a/set_mbtiles_metadata.py b/set_mbtiles_metadata.py new file mode 100644 index 0000000..723cf29 --- /dev/null +++ b/set_mbtiles_metadata.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 +import sqlite3 +import json +import subprocess +import sys +from pathlib import Path + +def read_gdal_bounds(raster): + """Legge i bounds WGS84 dal TIFF usando gdalinfo -json.""" + info = subprocess.check_output(["gdalinfo", "-json", raster]) + info = json.loads(info) + + if "wgs84Extent" in info: + coords = info["wgs84Extent"]["coordinates"][0] + minlon, minlat = coords[0] + maxlon, maxlat = coords[2] + else: + cc = info["cornerCoordinates"] + minlon, minlat = cc["lowerLeft"] + maxlon, maxlat = cc["upperRight"] + + # Normalizzazione latitudine + if minlat > maxlat: + minlat, maxlat = maxlat, minlat + + return (minlon, minlat, maxlon, maxlat) + + +def open_mbtiles(mbtiles): + """Apre MBTiles e garantisce che la tabella metadata esista.""" + conn = sqlite3.connect(mbtiles) + cur = conn.cursor() + + # Verifica tabella tiles + cur.execute("SELECT COUNT(*) FROM sqlite_master WHERE type='table' AND name='tiles';") + if cur.fetchone()[0] == 0: + raise RuntimeError("La tabella 'tiles' non esiste nell'MBTiles.") + + # Conta tile + cur.execute("SELECT COUNT(*) FROM tiles;") + tile_count = cur.fetchone()[0] + if tile_count == 0: + raise RuntimeError("La tabella 'tiles' Γ¨ vuota.") + + # Zoom min/max + cur.execute("SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;") + minzoom, maxzoom = cur.fetchone() + + # Crea metadata se manca + cur.execute(""" + CREATE TABLE IF NOT EXISTS metadata ( + name TEXT PRIMARY KEY, + value TEXT + ); + """) + + # Leggi metadata esistenti + cur.execute("SELECT name, value FROM metadata;") + metadata = dict(cur.fetchall()) + + return conn, cur, tile_count, minzoom, maxzoom, metadata + + +def write_metadata(cur, name, value): + cur.execute("REPLACE INTO metadata (name, value) VALUES (?, ?)", (name, value)) + + +def main(): + if len(sys.argv) < 3: + print("Uso: python3 set_mbtiles_metadata.py [name]") + sys.exit(1) + + raster = Path(sys.argv[1]) + mbtiles = Path(sys.argv[2]) + layer_name = sys.argv[3] if len(sys.argv) > 3 else mbtiles.stem + + if not raster.exists(): + print(f"Raster non trovato: {raster}") + sys.exit(1) + if not mbtiles.exists(): + print(f"MBTiles non trovato: {mbtiles}") + sys.exit(1) + + print("== Impostazione metadata MBTiles ==") + + # 1) Bounds dal TIFF + minlon, minlat, maxlon, maxlat = read_gdal_bounds(raster) + bounds = f"{minlon:.6f},{minlat:.6f},{maxlon:.6f},{maxlat:.6f}" + print(f"Bounds TIFF: {bounds}") + + # 2) Lettura MBTiles + conn, cur, tile_count, minzoom, maxzoom, metadata = open_mbtiles(mbtiles) + + print(f"Tile count: {tile_count}") + print(f"Zoom effettivi: {minzoom}..{maxzoom}") + + # 3) Center + ctr_lon = (minlon + maxlon) / 2 + ctr_lat = (minlat + maxlat) / 2 + ctr_zoom = (minzoom + maxzoom) // 2 + center = f"{ctr_lon:.6f},{ctr_lat:.6f},{ctr_zoom}" + + # 4) Scrittura metadata + write_metadata(cur, "name", layer_name) + write_metadata(cur, "format", "png") + write_metadata(cur, "type", "overlay") + write_metadata(cur, "scheme", "xyz") + write_metadata(cur, "bounds", bounds) + write_metadata(cur, "center", center) + write_metadata(cur, "minzoom", str(minzoom)) + write_metadata(cur, "maxzoom", str(maxzoom)) + + conn.commit() + conn.close() + + print("\n== Metadata scritti correttamente ==") + print(f"name: {layer_name}") + print(f"format: png") + print(f"type: overlay") + print(f"scheme: xyz") + print(f"bounds: {bounds}") + print(f"center: {center}") + print(f"minzoom: {minzoom}") + print(f"maxzoom: {maxzoom}") + print("\nMBTiles ora Γ¨ pronto per TileServer GL.") + + +if __name__ == "__main__": + main() + diff --git a/set_mbtiles_metadata_dir.py b/set_mbtiles_metadata_dir.py new file mode 100644 index 0000000..81f42ed --- /dev/null +++ b/set_mbtiles_metadata_dir.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +import sqlite3 +import json +import subprocess +import sys +from pathlib import Path + +def read_bounds(tif): + """Legge i bounds W,S,E,N dal TIFF usando gdalinfo -json, arrotondati a 6 decimali.""" + info = json.loads(subprocess.check_output(["gdalinfo", "-json", str(tif)])) + + if "wgs84Extent" in info: + coords = info["wgs84Extent"]["coordinates"][0] + minlon, minlat = coords[0] + maxlon, maxlat = coords[2] + else: + cc = info["cornerCoordinates"] + minlon, minlat = cc["lowerLeft"] + maxlon, maxlat = cc["upperRight"] + + # Normalizzazione + if minlon > maxlon: + minlon, maxlon = maxlon, minlon + if minlat > maxlat: + minlat, maxlat = maxlat, minlat + + return ( + round(minlon, 6), + round(minlat, 6), + round(maxlon, 6), + round(maxlat, 6), + ) + +def main(): + if len(sys.argv) != 3: + print("Uso: python3 set_mbtiles_metadata_multi.py ") + sys.exit(1) + + tif_dir = Path(sys.argv[1]) + mbtiles = Path(sys.argv[2]) + + if not tif_dir.exists() or not tif_dir.is_dir(): + print(f"Directory TIFF non valida: {tif_dir}") + sys.exit(1) + + if not mbtiles.exists(): + print(f"MBTiles non trovato: {mbtiles}") + sys.exit(1) + + # Trova tutti i TIFF nella directory + tifs = sorted([p for p in tif_dir.iterdir() if p.suffix.lower() in (".tif", ".tiff")]) + + if not tifs: + print("Nessun TIFF trovato nella directory.") + sys.exit(1) + + print(f"Trovati {len(tifs)} TIFF:") + for t in tifs: + print(" -", t.name) + + # Calcolo bounds unificati + bounds_list = [read_bounds(t) for t in tifs] + + W = min(b[0] for b in bounds_list) + S = min(b[1] for b in bounds_list) + E = max(b[2] for b in bounds_list) + N = max(b[3] for b in bounds_list) + + bounds = f"{W:.6f},{S:.6f},{E:.6f},{N:.6f}" + print("\nBounds unificati:", bounds) + + # Apri MBTiles + conn = sqlite3.connect(mbtiles) + cur = conn.cursor() + + cur.execute(""" + CREATE TABLE IF NOT EXISTS metadata ( + name TEXT PRIMARY KEY, + value TEXT + ); + """) + + # Zoom reali dal MBTiles + cur.execute("SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;") + minzoom, maxzoom = cur.fetchone() + + print(f"Zoom reali: {minzoom}..{maxzoom}") + + # Center + centerLon = (W + E) / 2 + centerLat = (S + N) / 2 + centerZoom = (minzoom + maxzoom) // 2 + center = f"{centerLon:.6f},{centerLat:.6f},{centerZoom}" + + # Scrittura metadata + cur.execute("REPLACE INTO metadata VALUES ('bounds', ?)", (bounds,)) + cur.execute("REPLACE INTO metadata VALUES ('center', ?)", (center,)) + cur.execute("REPLACE INTO metadata VALUES ('minzoom', ?)", (str(minzoom),)) + cur.execute("REPLACE INTO metadata VALUES ('maxzoom', ?)", (str(maxzoom),)) + cur.execute("REPLACE INTO metadata VALUES ('format', 'png')") + cur.execute("REPLACE INTO metadata VALUES ('type', 'overlay')") + cur.execute("REPLACE INTO metadata VALUES ('scheme', 'xyz')") + + conn.commit() + conn.close() + + print("\nMetadata scritti correttamente:") + print("bounds:", bounds) + print("center:", center) + print("minzoom:", minzoom) + print("maxzoom:", maxzoom) + print("format: png") + print("type: overlay") + print("scheme: xyz") + +if __name__ == "__main__": + main() + diff --git a/set_mbtiles_metadata_multi.py b/set_mbtiles_metadata_multi.py new file mode 100644 index 0000000..1dce4c3 --- /dev/null +++ b/set_mbtiles_metadata_multi.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +import sqlite3 +import json +import subprocess +import sys +from pathlib import Path + +def read_bounds(tif): + info = json.loads(subprocess.check_output(["gdalinfo", "-json", tif])) + if "wgs84Extent" in info: + coords = info["wgs84Extent"]["coordinates"][0] + minlon, minlat = coords[0] + maxlon, maxlat = coords[2] + else: + cc = info["cornerCoordinates"] + minlon, minlat = cc["lowerLeft"] + maxlon, maxlat = cc["upperRight"] + + if minlat > maxlat: + minlat, maxlat = maxlat, minlat + if minlon > maxlon: + minlon, maxlon = maxlon, minlon + + return ( + round(minlon, 6), + round(minlat, 6), + round(maxlon, 6), + round(maxlat, 6), + ) + +def main(): + if len(sys.argv) != 4: + print("Uso: python3 set_mbtiles_metadata_multi.py raster1.tif raster2.tif merged.mbtiles") + sys.exit(1) + + tif1 = Path(sys.argv[1]) + tif2 = Path(sys.argv[2]) + mbtiles = Path(sys.argv[3]) + + b1 = read_bounds(tif1) + b2 = read_bounds(tif2) + + W = min(b1[0], b2[0]) + S = min(b1[1], b2[1]) + E = max(b1[2], b2[2]) + N = max(b1[3], b2[3]) + + bounds = f"{W:.6f},{S:.6f},{E:.6f},{N:.6f}" + + conn = sqlite3.connect(mbtiles) + cur = conn.cursor() + + cur.execute("CREATE TABLE IF NOT EXISTS metadata (name TEXT PRIMARY KEY, value TEXT);") + + cur.execute("SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;") + minzoom, maxzoom = cur.fetchone() + + centerLon = (W + E) / 2 + centerLat = (S + N) / 2 + centerZoom = (minzoom + maxzoom) // 2 + center = f"{centerLon:.6f},{centerLat:.6f},{centerZoom}" + + cur.execute("REPLACE INTO metadata VALUES ('bounds', ?)", (bounds,)) + cur.execute("REPLACE INTO metadata VALUES ('center', ?)", (center,)) + cur.execute("REPLACE INTO metadata VALUES ('minzoom', ?)", (str(minzoom),)) + cur.execute("REPLACE INTO metadata VALUES ('maxzoom', ?)", (str(maxzoom),)) + cur.execute("REPLACE INTO metadata VALUES ('format', 'png')") + cur.execute("REPLACE INTO metadata VALUES ('type', 'overlay')") + cur.execute("REPLACE INTO metadata VALUES ('scheme', 'xyz')") + + conn.commit() + conn.close() + + print("Metadata scritti correttamente:") + print("bounds:", bounds) + print("center:", center) + print("minzoom:", minzoom) + print("maxzoom:", maxzoom) + print("format: png") + print("type: overlay") + print("scheme: xyz") + +if __name__ == "__main__": + main() + diff --git a/verify_and_fix_mbtiles.py b/verify_and_fix_mbtiles.py new file mode 100644 index 0000000..9ff896c --- /dev/null +++ b/verify_and_fix_mbtiles.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python3 +import sqlite3 +import json +import subprocess +import sys +from pathlib import Path + +def read_gdal_bounds(raster): + """Legge i bounds WGS84 dal TIFF usando gdalinfo -json.""" + info = subprocess.check_output(["gdalinfo", "-json", raster]) + info = json.loads(info) + + if "wgs84Extent" in info: + coords = info["wgs84Extent"]["coordinates"][0] + minlon, minlat = coords[0] + maxlon, maxlat = coords[2] + else: + cc = info["cornerCoordinates"] + minlon, minlat = cc["lowerLeft"] + maxlon, maxlat = cc["upperRight"] + + # Normalizzazione latitudine + if minlat > maxlat: + minlat, maxlat = maxlat, minlat + + return (minlon, minlat, maxlon, maxlat) + + +def read_mbtiles(mbtiles): + """Legge tile info e metadata da un MBTiles.""" + conn = sqlite3.connect(mbtiles) + cur = conn.cursor() + + # Verifica tabella tiles + cur.execute("SELECT COUNT(*) FROM sqlite_master WHERE type='table' AND name='tiles';") + if cur.fetchone()[0] == 0: + raise RuntimeError("La tabella 'tiles' non esiste nell'MBTiles.") + + # Conta tile + cur.execute("SELECT COUNT(*) FROM tiles;") + tile_count = cur.fetchone()[0] + if tile_count == 0: + raise RuntimeError("La tabella 'tiles' Γ¨ vuota.") + + # Zoom min/max + cur.execute("SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;") + minzoom, maxzoom = cur.fetchone() + + # Metadata + cur.execute("CREATE TABLE IF NOT EXISTS metadata (name TEXT PRIMARY KEY, value TEXT);") + cur.execute("SELECT name, value FROM metadata;") + metadata = dict(cur.fetchall()) + + return conn, cur, tile_count, minzoom, maxzoom, metadata + + +def write_metadata(cur, name, value): + cur.execute("REPLACE INTO metadata (name, value) VALUES (?, ?)", (name, value)) + + +def main(): + if len(sys.argv) != 3: + print("Uso: python3 verify_and_fix_mbtiles.py ") + sys.exit(1) + + raster = Path(sys.argv[1]) + mbtiles = Path(sys.argv[2]) + + if not raster.exists(): + print(f"Raster non trovato: {raster}") + sys.exit(1) + if not mbtiles.exists(): + print(f"MBTiles non trovato: {mbtiles}") + sys.exit(1) + + print("== Verifica e correzione metadata MBTiles ==") + + # 1) Bounds dal TIFF + minlon, minlat, maxlon, maxlat = read_gdal_bounds(raster) + expected_bounds = f"{minlon:.6f},{minlat:.6f},{maxlon:.6f},{maxlat:.6f}" + print(f"Bounds TIFF: {expected_bounds}") + + # 2) Lettura MBTiles + conn, cur, tile_count, minzoom, maxzoom, metadata = read_mbtiles(mbtiles) + + print(f"Tile count: {tile_count}") + print(f"Zoom effettivi: {minzoom}..{maxzoom}") + + # 3) Verifica/correzione bounds + mb_bounds = metadata.get("bounds") + if mb_bounds != expected_bounds: + print(f"✘ Bounds NON corretti β†’ FIX") + write_metadata(cur, "bounds", expected_bounds) + else: + print("βœ” Bounds OK") + + # 4) Verifica/correzione minzoom/maxzoom + mz = int(metadata.get("minzoom", -1)) + xz = int(metadata.get("maxzoom", -1)) + + if mz != minzoom or xz != maxzoom: + print(f"✘ minzoom/maxzoom NON coerenti β†’ FIX") + write_metadata(cur, "minzoom", str(minzoom)) + write_metadata(cur, "maxzoom", str(maxzoom)) + else: + print("βœ” minzoom/maxzoom OK") + + # 5) Verifica/correzione center + ctr_lon = (minlon + maxlon) / 2 + ctr_lat = (minlat + maxlat) / 2 + ctr_zoom = (minzoom + maxzoom) // 2 + expected_center = f"{ctr_lon:.6f},{ctr_lat:.6f},{ctr_zoom}" + + if metadata.get("center") != expected_center: + print(f"✘ Center NON corretto β†’ FIX") + write_metadata(cur, "center", expected_center) + else: + print("βœ” Center OK") + + # 6) Verifica/correzione format + if metadata.get("format") != "png": + print("✘ Format NON corretto β†’ FIX") + write_metadata(cur, "format", "png") + else: + print("βœ” Format OK") + + # 7) Verifica/correzione scheme + if metadata.get("scheme") != "xyz": + print("✘ Scheme NON corretto β†’ FIX") + write_metadata(cur, "scheme", "xyz") + else: + print("βœ” Scheme OK") + + # 8) Commit + conn.commit() + conn.close() + + print("\n== Correzione completata ==") + print("MBTiles ora ha metadata coerenti e validi.") + + +if __name__ == "__main__": + main() + diff --git a/verify_bounds_tileservergl.py b/verify_bounds_tileservergl.py new file mode 100644 index 0000000..8a6c1c0 --- /dev/null +++ b/verify_bounds_tileservergl.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +import sqlite3 +import json +import subprocess +import sys +from pathlib import Path + +EPS = 1e-6 # tolleranza per confronti float + +def almost_equal(a, b): + return abs(a - b) < EPS + +def read_gdal_bounds(raster): + """Legge i bounds W,S,E,N dal TIFF usando gdalinfo -json, arrotondati a 6 decimali.""" + info = subprocess.check_output(["gdalinfo", "-json", raster]) + info = json.loads(info) + + if "wgs84Extent" in info: + coords = info["wgs84Extent"]["coordinates"][0] + minlon, minlat = coords[0] + maxlon, maxlat = coords[2] + else: + cc = info["cornerCoordinates"] + minlon, minlat = cc["lowerLeft"] + maxlon, maxlat = cc["upperRight"] + + # Normalizzazione + if minlon > maxlon: + minlon, maxlon = maxlon, minlon + if minlat > maxlat: + minlat, maxlat = maxlat, minlat + + # Arrotondamento a 6 decimali (standard MBTiles) + return ( + round(minlon, 6), + round(minlat, 6), + round(maxlon, 6), + round(maxlat, 6), + ) + +def read_mbtiles_bounds(cur): + cur.execute("SELECT value FROM metadata WHERE name='bounds';") + row = cur.fetchone() + if not row: + return None + try: + w, s, e, n = map(float, row[0].split(",")) + return (w, s, e, n) + except: + return None + +def write_bounds(cur, bounds): + w, s, e, n = bounds + value = f"{w:.6f},{s:.6f},{e:.6f},{n:.6f}" + cur.execute("REPLACE INTO metadata (name,value) VALUES ('bounds', ?)", (value,)) + return value + +def main(): + if len(sys.argv) != 3: + print("Uso: python3 verify_bounds_tileservergl.py ") + sys.exit(1) + + raster = Path(sys.argv[1]) + mbtiles = Path(sys.argv[2]) + + if not raster.exists(): + print(f"Raster non trovato: {raster}") + sys.exit(1) + if not mbtiles.exists(): + print(f"MBTiles non trovato: {mbtiles}") + sys.exit(1) + + print("== Verifica bounds per TileServer GL (W,S,E,N) ==") + + # Bounds dal TIFF + tiff_bounds = read_gdal_bounds(raster) + print(f"Bounds TIFF (W,S,E,N): {tiff_bounds}") + + # Apri MBTiles + conn = sqlite3.connect(mbtiles) + cur = conn.cursor() + + # Assicura tabella metadata + cur.execute(""" + CREATE TABLE IF NOT EXISTS metadata ( + name TEXT PRIMARY KEY, + value TEXT + ); + """) + + # Bounds MBTiles + mb_bounds = read_mbtiles_bounds(cur) + print(f"Bounds MBTiles: {mb_bounds}") + + # Se mancano β†’ scrivili + if mb_bounds is None: + print("✘ Bounds mancanti β†’ FIX") + new_value = write_bounds(cur, tiff_bounds) + print(f"βœ” Bounds scritti: {new_value}") + conn.commit() + conn.close() + print("\n== Verifica completata ==") + return + + # Confronto con tolleranza + w1, s1, e1, n1 = tiff_bounds + w2, s2, e2, n2 = mb_bounds + + if not (almost_equal(w1, w2) and almost_equal(s1, s2) and almost_equal(e1, e2) and almost_equal(n1, n2)): + print("✘ Bounds NON corretti β†’ FIX") + new_value = write_bounds(cur, tiff_bounds) + print(f"βœ” Bounds aggiornati a: {new_value}") + conn.commit() + else: + print("βœ” Bounds OK (nessuna correzione necessaria)") + + conn.close() + print("\n== Verifica completata ==") + +if __name__ == "__main__": + main() + diff --git a/verify_mbtiles_metadata.py b/verify_mbtiles_metadata.py new file mode 100644 index 0000000..e21918f --- /dev/null +++ b/verify_mbtiles_metadata.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +import sqlite3 +import json +import subprocess +import sys +from pathlib import Path + +def read_gdal_bounds(raster): + """Legge i bounds WGS84 dal TIFF usando gdalinfo -json.""" + info = subprocess.check_output(["gdalinfo", "-json", raster]) + info = json.loads(info) + + if "wgs84Extent" in info: + coords = info["wgs84Extent"]["coordinates"][0] + minlon, minlat = coords[0] + maxlon, maxlat = coords[2] + else: + cc = info["cornerCoordinates"] + minlon, minlat = cc["lowerLeft"] + maxlon, maxlat = cc["upperRight"] + + # Normalizzazione latitudine + if minlat > maxlat: + minlat, maxlat = maxlat, minlat + + return (minlon, minlat, maxlon, maxlat) + + +def read_mbtiles_metadata(mbtiles): + """Legge metadata e tile info da un MBTiles.""" + conn = sqlite3.connect(mbtiles) + cur = conn.cursor() + + # Verifica tabella tiles + cur.execute("SELECT COUNT(*) FROM sqlite_master WHERE type='table' AND name='tiles';") + if cur.fetchone()[0] == 0: + raise RuntimeError("La tabella 'tiles' non esiste nell'MBTiles.") + + # Conta tile + cur.execute("SELECT COUNT(*) FROM tiles;") + tile_count = cur.fetchone()[0] + if tile_count == 0: + raise RuntimeError("La tabella 'tiles' Γ¨ vuota.") + + # Zoom min/max + cur.execute("SELECT MIN(zoom_level), MAX(zoom_level) FROM tiles;") + minzoom, maxzoom = cur.fetchone() + + # Metadata + cur.execute("SELECT name, value FROM metadata;") + metadata = dict(cur.fetchall()) + + conn.close() + return tile_count, minzoom, maxzoom, metadata + + +def main(): + if len(sys.argv) != 3: + print("Uso: python3 verify_mbtiles_metadata.py ") + sys.exit(1) + + raster = Path(sys.argv[1]) + mbtiles = Path(sys.argv[2]) + + if not raster.exists(): + print(f"Raster non trovato: {raster}") + sys.exit(1) + if not mbtiles.exists(): + print(f"MBTiles non trovato: {mbtiles}") + sys.exit(1) + + print("== Verifica metadata MBTiles ==") + + # 1) Bounds dal TIFF + minlon, minlat, maxlon, maxlat = read_gdal_bounds(raster) + print(f"Bounds TIFF: {minlon},{minlat},{maxlon},{maxlat}") + + # 2) Lettura MBTiles + tile_count, minzoom, maxzoom, metadata = read_mbtiles_metadata(mbtiles) + + print(f"Tile count: {tile_count}") + print(f"Zoom effettivi: {minzoom}..{maxzoom}") + + # 3) Confronto bounds + mb_bounds = metadata.get("bounds", None) + if mb_bounds: + mb_minlon, mb_minlat, mb_maxlon, mb_maxlat = map(float, mb_bounds.split(",")) + print(f"Bounds MBTiles: {mb_bounds}") + + if abs(mb_minlon - minlon) < 1e-4 and \ + abs(mb_minlat - minlat) < 1e-4 and \ + abs(mb_maxlon - maxlon) < 1e-4 and \ + abs(mb_maxlat - maxlat) < 1e-4: + print("βœ” Bounds OK") + else: + print("✘ Bounds NON corrispondono al TIFF") + else: + print("✘ Metadata 'bounds' mancante") + + # 4) Confronto minzoom/maxzoom + if "minzoom" in metadata and "maxzoom" in metadata: + mz = int(metadata["minzoom"]) + xz = int(metadata["maxzoom"]) + if mz == minzoom and xz == maxzoom: + print("βœ” minzoom/maxzoom OK") + else: + print(f"✘ minzoom/maxzoom NON coerenti: metadata={mz}..{xz}, reali={minzoom}..{maxzoom}") + else: + print("✘ Metadata minzoom/maxzoom mancanti") + + # 5) Center + if "center" in metadata: + print(f"Center: {metadata['center']}") + else: + print("✘ Metadata 'center' mancante") + + print("\n== Verifica completata ==") + + +if __name__ == "__main__": + main() + diff --git a/w.py b/w.py new file mode 100644 index 0000000..ecf7a43 --- /dev/null +++ b/w.py @@ -0,0 +1,68 @@ +import boto3 +import os +import logging +from botocore import UNSIGNED +from botocore.config import Config + +BUCKET = "copernicus-dem-30m" +OUTPUT_DIR = "./glo30_world" +LOG_FILE = "download_world.log" + +logging.basicConfig( + filename=LOG_FILE, + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", +) + +console = logging.getLogger("console") +console.setLevel(logging.INFO) +console.addHandler(logging.StreamHandler()) + +def log(msg): + logging.info(msg) + console.info(msg) + +os.makedirs(OUTPUT_DIR, exist_ok=True) + +s3 = boto3.client("s3", config=Config(signature_version=UNSIGNED)) + +log("Inizio scansione di tutte le cartelle globali...") + +folders = [] + +paginator = s3.get_paginator("list_objects_v2") + +# Legge SOLO le directory top-level +for page in paginator.paginate(Bucket=BUCKET, Delimiter="/"): + for prefix in page.get("CommonPrefixes", []): + folder = prefix["Prefix"] + + # Filtra solo i tile DEM + if folder.startswith("Copernicus_DSM_COG_10_") and +folder.endswith("_DEM/"): + folders.append(folder) + log(f"[FOUND] {folder}") + +log(f"Trovate {len(folders)} cartelle globali") + +# ----------------------------- +# DOWNLOAD SOLO IL DEM PRINCIPALE +# ----------------------------- +for folder in folders: + tif_name = folder[:-1] + ".tif" # rimuove "/" e aggiunge .tif + key = folder + tif_name.split("/")[-1] + out_path = os.path.join(OUTPUT_DIR, tif_name.split("/")[-1]) + + if os.path.exists(out_path): + log(f"[SKIP] {out_path} giΓ  presente") + continue + + log(f"[DOWNLOAD] {key}") + try: + s3.download_file(BUCKET, key, out_path) + log(f"[OK] {key} scaricato") + except Exception as e: + log(f"[ERROR] {key} β†’ {e}") + +log("Completato.") +