import laspy import datetime from collections import defaultdict from laspy.vlrs.known import WktCoordinateSystemVlr import os import folium from pyproj import Transformer # Dictionnaire des classes standard ASPRS class_names = { 1: "Non classé", 2: "Sol", 3: "Végétation basse", 4: "Végétation moyenne", 5: "Végétation haute", 6: "Bâtiment", 9: "Eau", 17: "Pont", 64: "Sursol pérenne", 65: "Artefact", 66: "Points virtuels (modélisation)", } # Demander à l'utilisateur de saisir le chemin et le nom du fichier file_path = input("Veuillez entrer le chemin complet et le nom du fichier .laz à analyser : ") file_name = os.path.basename(file_path) # Ouvrir le fichier .laz et afficher les informations à son sujet with laspy.open(file_path, laz_backend=laspy.LazBackend.Lazrs) as las_file: header = las_file.header point_format = header.point_format num_points = header.point_count point_attributes = point_format.dimension_names file_name_info = f"Nom du fichier : {file_name}" version_info = f"Version : {header.version}" num_points_info = f"Nombre de points : {num_points}" point_format_info = f"ID du format de point : {point_format.id}" point_attributes_info = "Attributs des points :
" + "
".join([f"- {attribute}" for attribute in point_attributes]) scales_info = f"Facteurs d'échelle : {header.scales}" offsets_info = f"Offsets : {header.offsets}" limits_info = (f"Limites :
" f" X : [{header.mins[0]}, {header.maxs[0]}]
" f" Y : [{header.mins[1]}, {header.maxs[1]}]
" f" Z : [{header.mins[2]}, {header.maxs[2]}]") # Examiner les VLR pour des informations de projection wkt_projection = None projection_info = "Projection : Non disponible" projection_info_full = "" vlrs = header.vlrs for vlr in vlrs: if isinstance(vlr, WktCoordinateSystemVlr): projection_info_full = vlr.string if len(projection_info_full) > 50: # Tronquer si trop long projection_info = "Projection : " + projection_info_full[:50] + "..." else: projection_info = "Projection : " + projection_info_full wkt_projection = vlr.string elif vlr.description.lower().find("projection") != -1 or vlr.description.lower().find("coordinate") != -1: projection_info_full = vlr.description if len(projection_info_full) > 50: # Tronquer si trop long projection_info = "Projection : " + projection_info_full[:50] + "..." else: projection_info = "Projection : " + projection_info_full # Lire les points points = las_file.read() gps_times = points.gps_time min_gps_time = min(gps_times) max_gps_time = max(gps_times) # Référence pour le temps GPS gps_epoch = datetime.datetime(2011, 9, 14, 0, 0, 0) start_time = gps_epoch + datetime.timedelta(seconds=min_gps_time) end_time = gps_epoch + datetime.timedelta(seconds=max_gps_time) start_time_info = f"Heure de début d'acquisition : {start_time}" end_time_info = f"Heure de fin d'acquisition : {end_time}" x_coords = points.x y_coords = points.y # Calculer les coins de la dalle x_min, x_max = min(x_coords), max(x_coords) y_min, y_max = min(y_coords), max(y_coords) # Calculer les dimensions de la dalle width = header.maxs[0] - header.mins[0] height = header.maxs[1] - header.mins[1] altitude_range = header.maxs[2] - header.mins[2] # Calculer la surface couverte surface_area = width * height # Calculer la densité de points point_density = num_points / surface_area # Trouver les points les plus bas et les plus hauts min_z = min(points.z) max_z = max(points.z) width_info = f"Largeur de la dalle (m) : {width}" height_info = f"Hauteur de la dalle (m) : {height}" altitude_range_info = f"Plage d'altitude de la dalle (m) : {altitude_range}" surface_area_info = f"Surface de la dalle (m²) : {surface_area}" point_density_info = f"Densité de points (points/m²) : {point_density}" min_z_info = f"Point le plus bas (Z) : {min_z}" max_z_info = f"Point le plus haut (Z) : {max_z}" # Compter les points par classe classifications = points.classification class_counts = defaultdict(int) for classification in classifications: class_counts[classification] += 1 class_info = "Classes et nombre de points associés :
" for class_number in sorted(class_counts): class_name = class_names.get(class_number, "User-defined or Reserved") count = class_counts[class_number] class_info += f" Classe {class_number} ({class_name}): {count} points
" # Afficher les informations print(file_name_info.replace("", "").replace("", "")) print(version_info.replace("", "").replace("", "")) print(num_points_info.replace("", "").replace("", "")) print(point_format_info.replace("", "").replace("", "")) print(point_attributes_info.replace("", "").replace("", "").replace("
", "\n")) print(scales_info.replace("", "").replace("", "")) print(offsets_info.replace("", "").replace("", "")) print(limits_info.replace("", "").replace("", "").replace("
", "\n")) print(projection_info_full) print(start_time_info.replace("", "").replace("", "")) print(end_time_info.replace("", "").replace("", "")) print(width_info.replace("", "").replace("", "")) print(height_info.replace("", "").replace("", "")) print(altitude_range_info.replace("", "").replace("", "")) print(surface_area_info.replace("", "").replace("", "")) print(point_density_info.replace("", "").replace("", "")) print(min_z_info.replace("", "").replace("", "")) print(max_z_info.replace("", "").replace("", "")) print(class_info.replace("", "").replace("", "").replace("
", "\n")) # Transformation des coordonnées des coins de la dalle if wkt_projection: transformer = Transformer.from_crs(wkt_projection, "EPSG:4326", always_xy=True) bottom_left = transformer.transform(x_min, y_min) bottom_right = transformer.transform(x_max, y_min) top_left = transformer.transform(x_min, y_max) top_right = transformer.transform(x_max, y_max) # Créer la carte avec folium m = folium.Map(location=[(bottom_left[1] + top_right[1]) / 2, (bottom_left[0] + top_right[0]) / 2], zoom_start=13) # Informations à afficher dans l'infobulle popup_content = ( file_name_info + "
" + version_info + "
" + num_points_info + "
" + point_format_info + "
" + point_attributes_info + "
" + scales_info + "
" + offsets_info + "
" + limits_info + "
" + projection_info + "
" + start_time_info + "
" + end_time_info + "
" + width_info + "
" + height_info + "
" + altitude_range_info + "
" + surface_area_info + "
" + point_density_info + "
" + min_z_info + "
" + max_z_info + "
" + class_info ) # Ajouter les coins de la dalle à la carte folium.Polygon( locations=[bottom_left[::-1], bottom_right[::-1], top_right[::-1], top_left[::-1], bottom_left[::-1]], color='blue', weight=2.5, fill=True, fill_color='blue', fill_opacity=0.2, popup=folium.Popup(popup_content, max_width=300) ).add_to(m) # Déterminer le répertoire de sortie pour le fichier HTML output_dir = os.path.dirname(file_path) output_path = os.path.join(output_dir, 'lidar_extent_map.html') # Enregistrer la carte m.save(output_path) print(f"Carte interactive créée : {output_path}")