relief/IGN.py

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])})"
	)