Code Python pour créer les Moyades en 3D
# -*- coding: utf-8 -*-
"""
Visualisation interactive du MNT autour de l’îlot des Moyades
Auteur : fcart
Date : juillet 2025
"""
# === IMPORTS ===
import pyvista as pv
import numpy as np
import rasterio
# === LECTURE DU FICHIER ===
path = r"C:\Users\bureau\geotif\0890_6235_Moyades\0890_6235\LITTO3D_FRA_0892_6234_20150529_LAMB93_RGF93_IGN69\MNT1m\LITTO3D_FRA_0892_6234_MNT_20150128_LAMB93_RGF93_IGN69.asc"
with rasterio.open(path) as src:
data = src.read(1)
transform = src.transform
nodata = src.nodata
# === MASQUAGE DES DONNÉES INEXISTANTES ===
data = np.where(data == nodata, np.nan, data)
# === CONSTRUCTION DES GRILLES X ET Y ===
nrows, ncols = data.shape
x0 = transform.c
y0 = transform.f
dx = transform.a
dy = transform.e # attention : négatif
x = np.arange(ncols) * dx + x0
y = np.arange(nrows) * dy + y0
xx, yy = np.meshgrid(x, y)
# === DÉFINITION MANUELLE DE LA ZONE D'INTÉRÊT (coordonnées Lambert 93) ===
# Ajustage des bornes selon besoin
xmin, xmax = 892550, 892950
#ymin, ymax = 6244002, 6246000
ymin, ymax = 6233600, 6233920
# === AIDE VISUELLE : APERÇU 2D AVEC MATPLOTLIB ===
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 8))
plt.imshow(data, extent=[xx.min(), xx.max(), yy.min(), yy.max()],
origin="lower", cmap="terrain", aspect='equal')
plt.colorbar(label="Altitude / Profondeur (m)")
plt.title("Aperçu de la zone MNT fusionnée (Lambert 93)")
plt.xlabel("X (Lambert 93)")
plt.ylabel("Y (Lambert 93)")
# Affichage de la zone d’intérêt actuelle (en rouge)
plt.axvline(xmin, color='red', linestyle='--', label='xmin / xmax')
plt.axvline(xmax, color='red', linestyle='--')
plt.axhline(ymin, color='red', linestyle='--', label='ymin / ymax')
plt.axhline(ymax, color='red', linestyle='--')
plt.legend(loc='lower right')
plt.grid(True)
plt.tight_layout()
plt.show()
# Création d'un masque logique pour découper la zone
mask_roi = (xx >= xmin) & (xx <= xmax) & (yy >= ymin) & (yy <= ymax)
row_indices, col_indices = np.where(mask_roi)
# Vérification de la présence de données valides dans la zone
if row_indices.size == 0 or col_indices.size == 0:
raise ValueError("⚠️ Aucune donnée trouvée dans la zone définie. Vérifie les bornes xmin/xmax/ymin/ymax.")
# Limites de découpe basées sur les indices valides
row_min, row_max = row_indices.min(), row_indices.max()
col_min, col_max = col_indices.min(), col_indices.max()
# === DÉCOUPAGE DES DONNÉES À LA ZONE D'INTÉRÊT ===
data_cropped = data[row_min:row_max+1, col_min:col_max+1]
xx_cropped = xx[row_min:row_max+1, col_min:col_max+1]
yy_cropped = yy[row_min:row_max+1, col_min:col_max+1]
# === STRUCTURED GRID PYVISTA ===
grid = pv.StructuredGrid(xx_cropped, yy_cropped, data_cropped)
grid["Profondeur"] = data_cropped.ravel(order="F")
# === SÉPARATION MER / TERRE ===
masque_mer = data_cropped <= 0
masque_terre = data_cropped > 0
data_mer = np.where(masque_mer, data_cropped, np.nan)
data_terre = np.where(masque_terre, data_cropped, np.nan)
grid_mer = pv.StructuredGrid(xx_cropped, yy_cropped, data_mer)
grid_mer["Profondeur"] = data_mer.ravel(order="F")
grid_terre = pv.StructuredGrid(xx_cropped, yy_cropped, data_terre)
grid_terre["Altitude"] = data_terre.ravel(order="F")
# === CONTOURS DE PROFONDEUR ===
niveaux = np.arange(-50, 1, 5)
contours = grid.contour(isosurfaces=niveaux, scalars="Profondeur")
# === CENTRAGE DE LA VUE ===
x_center = (xx_cropped.min() + xx_cropped.max()) / 2
y_center = (yy_cropped.min() + yy_cropped.max()) / 2
z_center = np.nanmax(data_cropped)
focus = (x_center, y_center, z_center)
camera_position = (x_center + 500, y_center + 500, z_center + 200)
viewup = (0, 0, 1)
# === AFFICHAGE INTERACTIF AVEC PYVISTA ===
plotter = pv.Plotter()
plotter.set_background("white") # ou "black" si tu préfères
plotter.add_mesh(grid_mer, cmap='viridis', scalars="Profondeur", nan_color="white", opacity=1.0)
plotter.add_mesh(grid_terre, color="saddlebrown", opacity=1.0)
plotter.add_mesh(contours, color="black", line_width=1)
#•plotter.add_scalar_bar(title="Profondeur (m)")
'''
plotter.add_mesh(parcours_line, color="cyan", line_width=3, label="Parcours de plongée")
plotter.add_mesh(bateau, color="red", label="Bateau")
plotter.add_legend()
'''
# === FLÈCHE DU NORD SIMPLIFIÉE ===
# Coordonnées dans le coin supérieur gauche du MNT
x_nord = xmin + 20
y_nord = ymax - 20
z_nord = -25 # Niveau de la mer
# Flèche simple : une ligne et une tête
nord_points = np.array([
[x_nord, y_nord, z_nord],
[x_nord, y_nord + 40, z_nord], # Tige vers le nord (Y+)
[x_nord - 5, y_nord + 35, z_nord], # Aile gauche
[x_nord + 5, y_nord + 35, z_nord], # Aile droite
[x_nord, y_nord + 40, z_nord], # Retour à la pointe
])
nord_line = pv.lines_from_points(nord_points[:2])
nord_head = pv.lines_from_points(np.vstack([nord_points[1], nord_points[2]])) + \
pv.lines_from_points(np.vstack([nord_points[1], nord_points[3]]))
# Ajout de la flèche au plotter
plotter.add_mesh(nord_line, color="black", line_width=3)
plotter.add_mesh(nord_head, color="black", line_width=3)
# Lettre "N"
plotter.add_point_labels(
[(x_nord, y_nord + 45, z_nord)],
["N"],
font_size=12,
point_color="white",
text_color="black",
shape_opacity=0.0,
always_visible=True
)
'''
# Étiquettes de profondeur (courbes de niveau principales)
for value in np.arange(-40, 1, 5):
try:
iso = contours.threshold((value - 0.1, value + 0.1), scalars="Profondeur")
if iso.n_points > 0:
point = iso.points[iso.n_points // 2]
plotter.add_point_labels(
[point],
[f"{int(value)} m"],
font_size=10,
point_color="white",
text_color="black",
shape_opacity=0.5,
always_visible=True
)
except Exception as e:
print(f"⚠️ Erreur pour la courbe {value} m : {e}")
'''
# Caméra
plotter.camera_position = (camera_position, focus, viewup)
plotter.camera.zoom(1.3)
# Affichage
plotter.show()
Commentaires
Enregistrer un commentaire