289 lines
12 KiB
Python
289 lines
12 KiB
Python
"""Téléchargement des dalles LiDAR HD de l'IGN pour les tuiles non générées.
|
||
|
||
Catalogue STAC (à jour) : https://browser.stac.teledetection.fr/collections/lidarhd
|
||
API : https://api.stac.teledetection.fr/collections/lidarhd/items
|
||
Fichiers (géoplateforme): https://data.geopf.fr/telechargement/download/...
|
||
|
||
Chaque dalle couvre 1 km × 1 km en Lambert 93 et est nommée par son coin
|
||
nord-ouest : LHD_FXX_{col}_{row}_PTS_LAMB93_IGN69.copc.laz
|
||
(col = X ouest en km, row = Y nord en km, cf. propriété STAC
|
||
"lidarhd:coordonnees_NW" au format "0816-6847").
|
||
"""
|
||
|
||
import json
|
||
import logging
|
||
import time
|
||
import urllib.parse
|
||
import urllib.request
|
||
from pathlib import Path
|
||
|
||
logger = logging.getLogger("lidar")
|
||
|
||
_STAC_ITEMS_URL = "https://api.stac.teledetection.fr/collections/lidarhd/items"
|
||
_HEADERS = {"User-Agent": "Mozilla/5.0 (lidar-archeo-pipeline)"}
|
||
|
||
|
||
def parse_tile_specs(args):
|
||
"""Convertit des spécifications "col,row" ou "col:row" en liste de tuples.
|
||
|
||
Args:
|
||
args: liste de chaînes (ex: ["1055,6882", "1056:6883"]).
|
||
|
||
Returns:
|
||
Liste de tuples (col, row).
|
||
|
||
Raises:
|
||
ValueError: si une spécification est mal formée.
|
||
"""
|
||
specs = []
|
||
for raw in args:
|
||
text = raw.strip().replace(":", ",").replace(";", ",")
|
||
parts = [p.strip() for p in text.split(",") if p.strip()]
|
||
if len(parts) != 2:
|
||
raise ValueError(f"Spécification de tuile invalide: {raw!r} (attendu: col,row)")
|
||
try:
|
||
col, row = int(parts[0]), int(parts[1])
|
||
except ValueError:
|
||
raise ValueError(f"Spécification de tuile invalide: {raw!r} (col et row doivent être des entiers)")
|
||
specs.append((col, row))
|
||
return specs
|
||
|
||
|
||
def tile_filename(col, row):
|
||
"""Nom de fichier LAZ standard d'une dalle (col, row)."""
|
||
return f"LHD_FXX_{col:04d}_{row:04d}_PTS_LAMB93_IGN69.copc.laz"
|
||
|
||
|
||
def _report_download(output_dir, tile, state, detail=None):
|
||
"""Émet un événement de téléchargement pour la file de génération (best-effort)."""
|
||
if output_dir is None:
|
||
return
|
||
from .progress import report_event
|
||
report_event(output_dir, tile, "download", state, detail=detail)
|
||
|
||
|
||
def _bbox_wgs84(col, row):
|
||
"""Bbox WGS84 de la dalle (col,row) pour la requête STAC (peut être élargie)."""
|
||
try:
|
||
from rasterio.warp import transform as warp_transform
|
||
xs = [col * 1000, (col + 1) * 1000, col * 1000, (col + 1) * 1000]
|
||
ys = [(row - 1) * 1000] * 2 + [row * 1000] * 2
|
||
lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', xs, ys)
|
||
except Exception:
|
||
lons = None
|
||
if lons is None:
|
||
try: # image légère (carte) : pyproj sans rasterio
|
||
from pyproj import Transformer
|
||
tr = Transformer.from_crs("EPSG:2154", "EPSG:4326", always_xy=True)
|
||
xs = [col * 1000, (col + 1) * 1000, col * 1000, (col + 1) * 1000]
|
||
ys = [(row - 1) * 1000] * 2 + [row * 1000] * 2
|
||
lons, lats = tr.transform(xs, ys)
|
||
except Exception:
|
||
lons = None
|
||
if lons is None:
|
||
from .index import _approx_l93_to_wgs84
|
||
pts = [_approx_l93_to_wgs84(x, y)
|
||
for x in (col * 1000, (col + 1) * 1000)
|
||
for y in ((row - 1) * 1000, row * 1000)]
|
||
lons = [p[0] for p in pts]
|
||
lats = [p[1] for p in pts]
|
||
pad = 0.005 # ~500 m de marge pour éviter les erreurs d'arrondi aux bords
|
||
return (min(lons) - pad, min(lats) - pad, max(lons) + pad, max(lats) + pad)
|
||
|
||
|
||
def match_feature(features, col, row):
|
||
"""Retourne l'item STAC correspondant à la dalle (col,row), sinon None."""
|
||
want = f"{col:04d}-{row:04d}"
|
||
for feature in features:
|
||
props = feature.get("properties", {})
|
||
if props.get("lidarhd:coordonnees_NW") == want:
|
||
return feature
|
||
return None
|
||
|
||
|
||
def find_tile_url(col, row, timeout=20, max_pages=5):
|
||
"""Cherche l'URL de téléchargement de la dalle (col,row) dans le catalogue STAC.
|
||
|
||
Returns:
|
||
URL (str) ou None si la dalle n'est pas (encore) publiée par l'IGN.
|
||
"""
|
||
feature = find_tile_feature(col, row, timeout=timeout, max_pages=max_pages)
|
||
return feature.get("assets", {}).get("data", {}).get("href") if feature else None
|
||
|
||
|
||
def tile_ign_metadata(col, row, timeout=10):
|
||
"""Fiche IGN d'une dalle : acquisition, capteurs, édition, téléchargement.
|
||
|
||
Returns:
|
||
dict normalisé (found=False si la dalle n'est pas publiée).
|
||
"""
|
||
feature = find_tile_feature(col, row, timeout=timeout)
|
||
if not feature:
|
||
return {"found": False}
|
||
p = feature.get("properties", {})
|
||
self_link = next((link.get("href") for link in feature.get("links", [])
|
||
if link.get("rel") == "self"), None)
|
||
return {
|
||
"found": True,
|
||
"download_url": feature.get("assets", {}).get("data", {}).get("href"),
|
||
"acquisition_start": p.get("start_datetime") or p.get("lidarhd:date_debut_acquisition"),
|
||
"acquisition_end": p.get("end_datetime") or p.get("lidarhd:date_fin_acquisition"),
|
||
"sensors": p.get("lidarhd:capteur") or [],
|
||
"mission": p.get("lidarhd:code_mission"),
|
||
"acquisition_operator": p.get("lidarhd:moe_acquisition"),
|
||
"acquisition_owner": p.get("lidarhd:moa_acquisition"),
|
||
"edition_date": p.get("lidarhd:date_edition"),
|
||
"classification_process": p.get("lidarhd:procede_classement"),
|
||
"points": p.get("lidarhd:nombre_points") or p.get("pc:count"),
|
||
"altimetry": p.get("lidarhd:systeme_altimetrique"),
|
||
"stac_url": self_link,
|
||
}
|
||
|
||
|
||
def find_tile_feature(col, row, timeout=20, max_pages=5):
|
||
"""Item STAC de la dalle (col,row), ou None si non publiée.
|
||
|
||
Suit la pagination du catalogue (lien rel=next) : sans elle, une dalle
|
||
au-delà de la première page de résultats semblait « introuvable ».
|
||
"""
|
||
w, s, e, n = _bbox_wgs84(col, row)
|
||
query = urllib.parse.urlencode({"bbox": f"{w:.6f},{s:.6f},{e:.6f},{n:.6f}", "limit": 50})
|
||
url = f"{_STAC_ITEMS_URL}?{query}"
|
||
for _ in range(max_pages):
|
||
req = urllib.request.Request(url, headers=_HEADERS)
|
||
with urllib.request.urlopen(req, timeout=timeout) as response:
|
||
data = json.loads(response.read().decode("utf-8"))
|
||
feature = match_feature(data.get("features", []), col, row)
|
||
if feature:
|
||
return feature
|
||
nxt = next((link.get("href") for link in data.get("links", [])
|
||
if link.get("rel") == "next" and link.get("href")), None)
|
||
if not nxt:
|
||
return None
|
||
url = nxt
|
||
return None
|
||
|
||
|
||
def download_file(url, dest_path, timeout=120, chunk=1024 * 1024):
|
||
"""Télécharge url vers dest_path en streaming (écriture atomique).
|
||
|
||
Écrit dans un .part puis renomme : un run interrompu (SIGKILL, panne)
|
||
ne laisse jamais un LAZ partiel qu'un run suivant considérerait comme
|
||
déjà téléchargé. Retourne la taille en octets.
|
||
"""
|
||
import os
|
||
dest_path = Path(dest_path)
|
||
tmp = dest_path.with_name(dest_path.name + ".part")
|
||
req = urllib.request.Request(url, headers=_HEADERS)
|
||
t0 = time.time()
|
||
try:
|
||
with urllib.request.urlopen(req, timeout=timeout) as response, open(tmp, "wb") as out:
|
||
done = 0
|
||
while True:
|
||
block = response.read(chunk)
|
||
if not block:
|
||
break
|
||
out.write(block)
|
||
done += len(block)
|
||
elapsed = time.time() - t0
|
||
logger.info(f" {done / 1e6:.0f} Mo en {elapsed:.0f}s"
|
||
f" ({done / 1e6 / max(elapsed, 0.1):.1f} Mo/s)")
|
||
os.replace(tmp, dest_path)
|
||
except BaseException:
|
||
tmp.unlink(missing_ok=True)
|
||
raise
|
||
return done
|
||
|
||
|
||
def fetch_tiles(input_dir, specs, output_dir=None, only_viz=None, resolutions=(0.5,), force=False):
|
||
"""Télécharge les dalles IGN spécifiées, sauf celles inutiles au traitement.
|
||
|
||
Args:
|
||
input_dir: dossier des fichiers LAZ (écriture autorisée requise).
|
||
specs: liste de tuples (col, row).
|
||
output_dir: dossier de sortie (optionnel) — permet d'ignorer les
|
||
tuiles déjà complètes.
|
||
only_viz: noms d'étapes des visualisations demandées (ex: ['aspect']).
|
||
Si fourni, une tuile n'est ignorée que si elle possède déjà
|
||
toutes ces visualisations aux résolutions demandées —
|
||
sinon son LAZ est (re)téléchargé pour compléter les
|
||
visualisations manquantes.
|
||
resolutions: résolutions attendues (m/px) pour considérer une tuile
|
||
complète quand only_viz est fourni.
|
||
force: True pour télécharger même les tuiles complètes (régénération :
|
||
leur LAZ est requis pour retraiter avec --force).
|
||
|
||
Returns:
|
||
Liste des chemins téléchargés.
|
||
"""
|
||
from .dtm import EDGE_NEIGHBORS_DIRNAME
|
||
input_dir = Path(input_dir)
|
||
complete = None
|
||
if output_dir is not None and only_viz and not force:
|
||
from .index import cells_with_all_viz, step_to_keyword
|
||
keys = [step_to_keyword(v) for v in only_viz]
|
||
complete = cells_with_all_viz(Path(output_dir) / "visualisations",
|
||
keys, resolutions)
|
||
# Explicite : jamais de promotion quand on télécharge directement DANS
|
||
# edge_neighbors/ (appel de _fetch_edge_neighbors avec input_dir=edge_dir)
|
||
# — sinon une dalle voisine pourrait se retrouver déplacée vers elle-même.
|
||
is_edge_dir = input_dir.name == EDGE_NEIGHBORS_DIRNAME
|
||
downloaded = []
|
||
for col, row in specs:
|
||
name = tile_filename(col, row)
|
||
tile = name[:-len(".copc.laz")] if name.endswith(".copc.laz") else name
|
||
dest = input_dir / name
|
||
if dest.exists():
|
||
logger.info(f" {name} : déjà présent dans input/ — aucun téléchargement")
|
||
_report_download(output_dir, tile, "skip", "déjà dans input/")
|
||
continue
|
||
if not is_edge_dir:
|
||
# dest.exists() ci-dessus ne teste que le nom exact : un ".part"
|
||
# voisin n'y correspond jamais, donc .exists() suffit à exclure
|
||
# un téléchargement voisin encore en cours.
|
||
neighbor_copy = input_dir / EDGE_NEIGHBORS_DIRNAME / name
|
||
if neighbor_copy.exists():
|
||
import os
|
||
os.replace(neighbor_copy, dest)
|
||
logger.info(f" {name} : déjà téléchargée comme voisine "
|
||
f"— déplacée dans input/")
|
||
_report_download(output_dir, tile, "skip",
|
||
"déjà téléchargée comme voisine — déplacée dans input/")
|
||
continue
|
||
if output_dir is not None and not force:
|
||
if complete is not None:
|
||
if (col, row) in complete:
|
||
logger.info(f" {name} : visualisations déjà complètes — ignorée")
|
||
_report_download(output_dir, tile, "skip", "visualisations déjà complètes")
|
||
continue
|
||
else:
|
||
vis_dir = Path(output_dir) / "visualisations"
|
||
if list(vis_dir.glob(f"LHD_FXX_{col:04d}_{row:04d}_PTS*")):
|
||
logger.info(f" {name} : visualisations déjà générées — ignorée")
|
||
_report_download(output_dir, tile, "skip", "visualisations déjà générées")
|
||
continue
|
||
logger.info(f" {name} : recherche dans le catalogue IGN...")
|
||
_report_download(output_dir, tile, "start", "catalogue IGN")
|
||
try:
|
||
url = find_tile_url(col, row)
|
||
except Exception as e:
|
||
logger.warning(f" ✗ {name} : erreur catalogue ({e})")
|
||
_report_download(output_dir, tile, "fail", f"erreur catalogue : {e}")
|
||
continue
|
||
if not url:
|
||
logger.warning(f" ✗ {name} : introuvable dans le catalogue IGN (zone non publiée ?)")
|
||
_report_download(output_dir, tile, "fail", "introuvable dans le catalogue IGN")
|
||
continue
|
||
logger.info(f" {name} : téléchargement depuis la géoplateforme...")
|
||
_report_download(output_dir, tile, "start", "géoplateforme")
|
||
try:
|
||
download_file(url, dest)
|
||
logger.info(f" ✓ {name} téléchargée")
|
||
_report_download(output_dir, tile, "ok")
|
||
downloaded.append(dest)
|
||
except Exception as e:
|
||
dest.unlink(missing_ok=True)
|
||
logger.warning(f" ✗ {name} : échec du téléchargement ({e})")
|
||
_report_download(output_dir, tile, "fail", f"échec du téléchargement : {e}")
|
||
return downloaded
|