import os import json import tempfile import py7zr import numpy as np from scipy.ndimage import zoom def Construit(nom, listeLambert, dossier, reduction=1): """ Construit un fichier RCL contenant une matrice d'altitudes. Args: nom: Nom du fichier de sortie sans extension. listeLambert: Liste des tuiles Lambert classées Nord-Ouest vers Sud-Est. dossier: Dossier contenant les archives RGEALTI .asc.7z. reduction: Facteur de réduction de résolution. 1 = résolution originale 2 = réduction par 2 avec interpolation bicubique 5 = réduction par 5 avec interpolation bicubique """ if reduction < 1: raise ValueError("La réduction doit être >= 1") nom = nom + ".rcl" taille = int(len(listeLambert) ** 0.5) if taille * taille != len(listeLambert): raise ValueError("La liste ne forme pas un carré complet") resultat = [] with tempfile.TemporaryDirectory() as tmp: # Parcours des lignes de tuiles du Nord vers le Sud for ligne_tuile in range(taille): ligne_fichiers = listeLambert[ ligne_tuile * taille: (ligne_tuile + 1) * taille ] fichiers = [] # Ouverture des tuiles Est-Ouest for lambert in ligne_fichiers: archive_path = os.path.join( dossier, f"RGEALTI_FXX_{lambert.replace(';', '_')}_MNT_LAMB93_IGN69.asc.7z" ) with py7zr.SevenZipFile( archive_path, mode="r" ) as archive: noms = archive.getnames() archive.extractall(path=tmp) nom_asc = None for fichier in noms: if fichier.endswith(".asc"): nom_asc = os.path.join(tmp, fichier) break if nom_asc is None: raise FileNotFoundError( f"Aucun fichier ASC trouvé dans {archive_path}" ) f = open(nom_asc, "r") # Suppression de l'entête ASC for _ in range(6): next(f) fichiers.append(f) # Lecture des 1000 lignes de cette bande for _ in range(1000): ligne_complete = [] for f in fichiers: ligne_complete.extend( float(x) for x in f.readline().split() ) resultat.append(ligne_complete) for f in fichiers: f.close() # Conversion numpy terrain = np.array( resultat, dtype=np.float32 ) # Réduction bicubique if reduction > 1: nouvelle_hauteur = terrain.shape[0] / reduction nouvelle_largeur = terrain.shape[1] / reduction terrain = zoom( terrain, ( nouvelle_hauteur / terrain.shape[0], nouvelle_largeur / terrain.shape[1] ), order=3 ) # Retour en liste pour JSON resultat = terrain.tolist() with open(nom, "w") as f: json.dump(resultat, f) print( f"RCL créé : {nom} " f"({len(resultat)} x {len(resultat[0])})" )