Je suis ravi de partager un projet qui redéfinit l’extraction et la visualisation des données géospatiales depuis OpenStreetMap, ainsi que le traitement de données d’élévation et des images satellites. Ce script Python, conçu pour les professionnels et les amateurs de géodonnées, offre une suite complète de fonctions pour une analyse approfondie et une visualisation immersive.
Pour commencer, nous avons tout d’abord installer les bibliothèques et les dépendances nécessaires puis nous procéderons à l’extraction des données OSM et à l’attribution des couleurs…
[ ]
!pip install osmnx geopy networkx osmium folium geopandas matplotlib contextily ezdxf
[ ]
import ssl
import certifi
import geopy
import osmnx as ox
import networkx as nx
from PIL import Image
import folium
from pyproj import Proj, transform
import pandas as pd
import geopandas as gpd
import contextily as ctx
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from shapely.geometry import Point
import plotly.graph_objects as go
import numpy as np
import ezdxf
import rasterio
# Variables pour extraire les données de la map:
address = 'fès, maroc'
ctx = ssl.create_default_context(cafile=certifi.where())
geopy.geocoders.options.default_ssl_context = ctx
geolocator = geopy.geocoders.Nominatim(timeout=3, user_agent = 'test')
location = geolocator.geocode(address, addressdetails=True)
location = location.raw
center_point = (float(location['lat']),float(location['lon']))
latitude = float(location['lat'])
longitude = float(location['lon'])
dist = 3000
bbox = ox.utils_geo.bbox_from_point(center_point, dist=dist)
[ ]
tags = {"building" : True}
# Extraire les éléments requis (bâtiments, réseau routier, cours d'eau, forêts et espaces vert)
G_Roads = ox.graph_from_point(center_point, dist=dist, retain_all=True, simplify = True, network_type='all')
G_Water = ox.graph_from_point(center_point, dist=dist, dist_type='bbox', network_type='all', simplify=True, retain_all=True, clean_periphery=False, custom_filter=["waterway"])
G_Green = ox.graph_from_point(center_point, dist=dist, dist_type='bbox', network_type='all',
simplify=True, retain_all=True, truncate_by_edge=False,
clean_periphery=False, custom_filter=["natural"])
GDF = ox.features_from_point(center_point, dist=dist, tags=tags)
# Extraire les points d'intérêt
GPOIs = ox.features_from_point(center_point, dist=dist, tags={"amenity": True}) [ ]
def TextCoords(coords):
"""
Cette fonction prend des coordonnées géographiques (latitude et longitude) en entrée
et renvoie une chaîne de caractères formatée indiquant la direction cardinale (Nord, Sud, Est ou Ouest).
Args:
coords (tuple): Un tuple contenant les coordonnées (latitude, longitude).
Returns:
str: La chaîne de caractères formatée avec les directions cardinales.
"""
out = ''
# Formate la latitude
out = out + str(round(abs(coords[0]), 5)) + 'º '
if coords[0] > 0:
out = out + 'N' # Nord
else:
out = out + 'S' # Sud
out = out + ' ' # Espacement
# Formate la longitude
out = out + str(round(abs(coords[1]), 5)) + 'º '
if coords[1] > 0:
out = out + 'E' # Est
else:
out = out + 'W' # Ouest
return out
[ ]
# Maintenant nous pouvons commencer à préparer les données pour la visualisation:
figure_size = (5.5,5.5) # Pour générer une image carrée
fig_dpi = 600 # Choix de résolution pour la qualité d'image, 300 est la norme pour l'impression
coordinate_text = TextCoords(center_point)
title = 'FÈS'
subtitle = 'MAROC'
text_color = (237, 230, 211)
# Configurer les couleurs:
# Classique
width_scale= 0.5
road_color_set = ["#a6a6a6" ,"#676767" ,"#454545" ,"#bdbdbd" ,"#d5d5d5" ]
road_linewidth_set = [0.1/width_scale,0.15/width_scale,0.25/width_scale,0.35/width_scale,0.45/width_scale]
highway_color= "#ffff";highway_linewidth= 0.5/width_scale
bgcolor = "#061529"
water_color = '#72b1b1'
water_linewidth = 0.5/width_scale
green_color = '#005f73'
green_linewidth = 0.5/width_scale
footprint_color = 'orange'
point_color = '#e8e8e4'
point_size = 20
# # Technique
# road_color_set = ["#6c757d" ,"#adb5bd" ,"#ced4da" ,"#dee2e6" ,"#e9ecef" ]
# road_linewidth_set = [0.1/width_scale,0.15/width_scale,0.25/width_scale,0.35/width_scale,0.45/width_scale]
# highway_color= "#343a40";highway_linewidth= 0.5/width_scale
# bgcolor = "#fffffc"
# water_color = 'blue'
# water_linewidth = 0.5/width_scale
# green_color = 'green'
# green_linewidth = 0.5/width_scale
# footprint_color = 'purple'
# point_color = '#f5cb5c'
# point_size = 30
[ ]
# Générer le tracé du réseau routier
# Déballer les données pour appliquer les conditions de couleurs
u = [];v = [];key = [];data = []
for uu, vv, kkey, ddata in G_Roads.edges(keys=True, data=True):
u.append(uu);v.append(vv);key.append(kkey);data.append(ddata)
# Appliquer les couleurs au réseau routier
roadColors = [];roadWidths = []
for item in data:
if "length" in item.keys():
if item["length"] <= 50:linewidth = road_linewidth_set[0];color = road_color_set[0]
elif item["length"] > 50 and item["length"] <= 100:linewidth = road_linewidth_set[1];color = road_color_set[1]
elif item["length"] > 100 and item["length"] <= 150:linewidth = road_linewidth_set[2];color = road_color_set[2]
elif item["length"] > 150 and item["length"] <= 200:color = road_color_set[3];linewidth = road_linewidth_set[3]
else:
color = road_color_set[4];linewidth = road_linewidth_set[4]
if "primary" in item["highway"]:
linewidth = highway_linewidth;color = highway_color
else:
color = road_color_set[0];linewidth = road_linewidth_set[0]
roadColors.append(color)
roadWidths.append(linewidth)
À présent, nous sommes prêts pour commencer à traiter les données extraites.
Dans la partie précédente du code, nous avons procédé à l’extraction de 4 types de données. Les données du réseau routier, les données des cours d’eau, les données des espaces verts et enfin celles des bâtiments. Dans cette partie, nous allons tout simplement générer des images heute résolution pour chaque type de données, et puis nous allons combiner ces images pour obtenir une vue d’ensemble des données.
[ ]
# Enfin nous pouvons commencer à générer nos images:
fig, ax = ox.plot_graph(G_Roads, node_size=0,
dpi = fig_dpi,bgcolor = bgcolor, figsize=figure_size,bbox=bbox,
save = False,show=False ,edge_color=roadColors,
edge_linewidth=roadWidths, edge_alpha=1)
fig.tight_layout(pad=0)
fig.savefig("BASE.png", dpi=fig_dpi, format="png", bbox_inches='tight',
facecolor=bgcolor, transparent=False,edgecolor=bgcolor,pad_inches = 0)
[ ]
fig, ax = ox.plot_graph(G_Water, node_size=0,figsize=figure_size,bbox=bbox,
dpi = fig_dpi, save = False, show=False, edge_color=water_color,
edge_linewidth=water_linewidth, edge_alpha=1)
fig.tight_layout(pad=0)
fig.savefig("WATER.png", dpi=fig_dpi, format="png", bbox_inches='tight',
facecolor=bgcolor, transparent=True,edgecolor=bgcolor,pad_inches = 0)
[ ]
fig, ax = ox.plot_graph(G_Green, node_size=0,figsize=figure_size,bbox=bbox,
dpi = fig_dpi, save = False, show=False, edge_color=green_color,
edge_linewidth=green_linewidth, edge_alpha=1)
fig.tight_layout(pad=0)
fig.savefig("GREEN.png", dpi=fig_dpi, format="png", bbox_inches='tight',
facecolor=bgcolor, transparent=True,edgecolor=bgcolor,pad_inches = 0)
[ ]
fig, ax = ox.plot_footprints(GDF, ax=None, figsize=figure_size, color=footprint_color,bbox=bbox, bgcolor='#111111',save=False, show=False, close=False, filepath=None, dpi=fig_dpi)
# Ajout des POI
ax.scatter(center_point[1],center_point[0], c=point_color, s=point_size)
fig.tight_layout(pad=0)
fig.savefig("FOOTPRINT.png", dpi=fig_dpi, format="png", bbox_inches='tight',
facecolor=bgcolor, transparent=True,edgecolor=bgcolor,pad_inches = 0)
[ ]
# Accéder aux images PNG transparentes (remplacez les noms de fichier par les vôtres)
image1 = Image.open('BASE.png')
image2 = Image.open('WATER.png')
image3 = Image.open('GREEN.png')
image4 = Image.open('FOOTPRINT.png')
# Créer une nouvelle image vide avec la taille de l'une des images
largeur, hauteur = image1.size
image_combinee = Image.new('RGBA', (largeur, hauteur))
# Superposer les images (vous pouvez ajuster les coordonnées selon vos besoins)
image_combinee.paste(image1, (0, 0), image1)
image_combinee.paste(image2, (0, 0), image2)
image_combinee.paste(image3, (0, 0), image3)
image_combinee.paste(image4, (0, 0), image4)
# Enregistrer l'image composite
image_combinee.save('COMBI.png')
image_combinee.save("COMBI_georeferenced.tiff", format="TIFF")
display(image_combinee)
Et maintenant, passons aux choses plus sérieuses et créons une carte interactive pour chaque type de données. Le réseau routier, les cours d’eau, les espaces verts et enfin les bâtiments…
[ ]
# Créer une carte centrée sur le point d'intérêt
m = folium.Map(location=center_point, zoom_start=14)
# Ajouter des lignes pour le réseau routier
for u, v, data in G_Roads.edges(data=True):
start_point = G_Roads.nodes[u]['y'], G_Roads.nodes[u]['x']
end_point = G_Roads.nodes[v]['y'], G_Roads.nodes[v]['x']
folium.PolyLine([start_point, end_point], color='blue', weight=2, opacity=1).add_to(m)
# Sauvegarder la carte dans un fichier HTML
m.save('road_network_map.html')
m
[ ]
# Créer une carte centrée sur le point d'intérêt
m = folium.Map(location=center_point, zoom_start=14)
# Ajouter des lignes pour les cours d'eau
for u, v, data in G_Water.edges(data=True):
start_point = G_Water.nodes[u]['y'], G_Water.nodes[u]['x']
end_point = G_Water.nodes[v]['y'], G_Water.nodes[v]['x']
folium.PolyLine([start_point, end_point], color='blue', weight=2, opacity=1).add_to(m)
# Sauvegarder la carte dans un fichier HTML
m.save('waterway_network_map.html')
m
[ ]
# Créer une carte centrée sur le point d'intérêt
#center_point = [34.0383594, -5.0043808]
m = folium.Map(location=center_point, zoom_start=16)
# Ajouter des marqueurs pour les bâtiments avec les coordonnées du centre
popup_text = "<b>Espace Vert</b>"
for node in G_Green.nodes():
if 'y' in G_Green.nodes[node] and 'x' in G_Green.nodes[node]:
location = [G_Green.nodes[node]['y'], G_Green.nodes[node]['x']]
else:
continue # Skip nodes without 'lat' and 'lon' keys.
folium.Marker(
location=location,
popup=popup_text,
icon=folium.Icon(color='orange')
).add_to(m)
# Sauvegarder la carte dans un fichier HTML
m.save('map.html')
# Display the map
m
Et ainsi, pour la carte interactive des bâtiments, nous allons rajouter les coordonnées du centre de chaque bâtiment dans un pop-up qui apparaît en cliquant sur le marqueur du bâtiment, cela nous servira à extraire un bâtiment pour les étapes suivantes du traitement…
[ ]
# Créer une carte centrée sur le point d'intérêt
m = folium.Map(location=center_point, zoom_start=19)
# Ajouter des marqueurs pour les bâtiments avec les coordonnées du centre
for idx, row in GDF.iterrows():
# Calculer les coordonnées du centre du bâtiment
center = row.geometry.centroid
popup_text = f"Centre: {center.y:.5f}, {center.x:.5f}"
folium.Marker(
location=[center.y, center.x],
popup=popup_text,
icon=folium.Icon(color='orange')
).add_to(m)
# Sauvegarder la carte dans un fichier HTML
m.save('map.html')
m
Mais avant d’aller plus loin, définissons une petite fonction qui nous servira à transformer les projections pour effectuer des calculs de distance et de superficie. Nous utilisons la fonction pour calculer la densité par kilomètre carré dans la zone d’étude…
[ ]
# Définir la projection WGS84 (système de coordonnées géographiques)
wgs84 = Proj('epsg:4326')
# Définir la projection UTM pour la zone 30N
utm = Proj('epsg:26191')
# Fonction pour convertir des coordonnées de degrés en mètres
def convert_coords(lat, lon):
x, y = transform(wgs84, utm, lon, lat)
return x, y
[ ]
# Convertir le graphe en GeoDataFrames pour les nœuds et les arêtes
gdf_nodes, gdf_edges = ox.graph_to_gdfs(G_Roads)
# 'convert_coords' est notre fonction de conversion de coordonnées
gdf_nodes['x_m'], gdf_nodes['y_m'] = zip(*gdf_nodes.apply(lambda row: convert_coords(row['y'], row['x']), axis=1))
# Créer un polygone avec les nouvelles coordonnées en mètres
coords_m = [(data['x_m'], data['y_m']) for _, data in gdf_nodes.iterrows()]
polygon_m = Polygon(coords_m)
area_m = polygon_m.area
# Convertir la superficie en kilomètres carrés
area_sq_km = area_m / 1e6
# Notre fonction de densité
def calculate_density(gdf, area_sq_km):
return len(gdf) / area_sq_km
# Ensuite, vous pouvez calculer la densité comme avant
densities = {
'Bâtiments': calculate_density(GDF, area_sq_km), 'Routes': calculate_density(gdf_edges, area_sq_km),
'Cours d'eau': calculate_density(G_Water, area_sq_km),
'Espaces verts': calculate_density(G_Green, area_sq_km)
}
# Afficher les densités
print("Densités par kilomètre carré:")
for feature, density in densities.items():
print(f"{feature}: {density:.2f}")
# Fonction pour visualiser le réseau routier
def plot_graph(G):
fig, ax = plt.subplots(figsize=(12, 8))
nodes = nx.draw_networkx_nodes(G, pos=nx.spring_layout(G), node_color='black', node_size=5)
edges = nx.draw_networkx_edges(G, pos=nx.spring_layout(G), edge_color='gray')
plt.show()
# Fonction pour visualiser les GeoDataFrames
def plot_gdfs(gdfs, colors, labels):
fig, ax = plt.subplots(figsize=(12, 8))
for gdf, color, label in zip(gdfs, colors, labels):
gdf.plot(ax=ax, color=color, label=label)
plt.legend()
plt.show()
# Visualiser les données
#plot_graph(G_Roads)
#plot_gdfs([GDF, gdf_edges],
#['red', 'gray'],
#['Bâtiments', 'Routes'])
Densités par kilomètre carré: Bâtiments: 32.36 Routes: 246.79 Cours d'eau: 0.12 Espaces verts: 2.81
Et voici une image géoréférencée sur la nouvelle projection (EPSG26191), représentant tous les bâtiments dans la zone d’étude.
[ ]
# Créer une figure et un axe avec Matplotlib
fig, ax = plt.subplots(figsize=(5, 5))
# Dessiner les données sur l'axe
GDF.to_crs(epsg=26191).plot(ax=ax)
# Enregistrer la figure avec le fond de carte géoréférencé
plt.savefig('carte_geo.tiff', dpi=300)
Vous avez envie d’aller plus loin ! ? Dans ce cas, pourquoi ne pas essayer d’extraire un bâtiment spécifié par ses coordonnées, que nous pouvons obtenir depuis la carte interactive des bâtiments que nous avons créé précédemment, pour ainsi le visualiser et calculer sa superficie et son périmètre…
[ ]
# Assurez-vous que 'GDF' est votre GeoDataFrame contenant les données des bâtiments
# Reprojeter le GeoDataFrame en UTM
#GDF2 = GDF.to_crs(epsg=26191)
# Coordonnées du bâtiment à visualiser avec une tolérance
latitude_batiment = 34.03448 # Remplacez par la latitude du bâtiment
longitude_batiment = -5.01622 # Remplacez par la longitude du bâtiment
tolerance = 0.0001 # Tolérance spatiale, par exemple 10 mètres
# Créer un point avec les coordonnées du bâtiment
point_batiment = Point(longitude_batiment, latitude_batiment)
# Sélectionner le bâtiment dans le GeoDataFrame avec une tolérance
batiment = GDF[GDF.distance(point_batiment) <= tolerance]
batiment2 = batiment.to_crs(epsg=26191)
# Vérifier si le bâtiment est trouvé
if not batiment2.empty:
# Créer une figure pour la visualisation
fig, ax = plt.subplots()
# Tracer le contour du bâtiment
batiment2.plot(ax=ax, color='blue')
# Calculer et afficher la superficie et la longueur des côtés
superficie = batiment2.geometry.area.iloc[0] # Superficie en mètres carrés
perimetre = batiment2.geometry.length.iloc[0] # Périmètre en mètres
# Afficher les informations sur l'image
info_text = f"Superficie: {superficie:.2f} m²nPérimètre: {perimetre:.2f} m"
plt.text(0.05, 0.95, info_text, transform=ax.transAxes, fontsize=12,
verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.5))
# Afficher l'image
plt.show()
else:
print("Aucun bâtiment trouvé à proximité de ces coordonnées.")
Nous allons ensuite visualiser le bâtiment extrait dans une forme tridimensionnelle statique représentée, dans un premier temps, par les coordonnées extérieures du bâtiment, puis nous allons le visualiser dans une forme représentée par tous les points de son polygone…
[ ]
# Vérifier si le bâtiment est trouvé
if not batiment2.empty:
# Créer une figure 3D pour la visualisation
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Obtenir les coordonnées x, y du polygone
x, y = batiment2.geometry.iloc[0].exterior.xy
# Spécifier une hauteur pour le bâtiment
hauteur = 50 # Hauteur en mètres, remplacez par la hauteur réelle si disponible
# Créer un tableau de la même hauteur pour chaque point
z = np.full(len(x), hauteur)
# Tracer le contour du bâtiment en 3D
ax.plot_trisurf(x, y, z, color='blue', alpha=0.7)
# Afficher l'image
plt.show()
else:
print("Aucun bâtiment trouvé à proximité de ces coordonnées.")
[ ]
# Vérifier si le bâtiment est trouvé
if not batiment2.empty:
# Créer une figure 3D pour la visualisation
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Obtenir les coordonnées x, y du polygone
coords = batiment2.geometry.iloc[0].exterior.coords
x, y = zip(*coords)
# Spécifier une hauteur pour le bâtiment
hauteur = 50 # Hauteur en mètres, remplacez par la hauteur réelle si disponible
# Créer les points supérieurs du bâtiment en ajoutant la hauteur
z = [hauteur] * len(x)
# Combinez les coordonnées inférieures et supérieures
verts = [list(zip(x, y, [0]*len(x))), list(zip(x, y, z))]
# Créer les faces du polygone 3D
faces = Poly3DCollection(verts, alpha=0.7, facecolors='cyan', linewidths=1, edgecolors='r', zsort='min')
# Ajouter les faces au graphique
ax.add_collection3d(faces)
# Définir les limites de la figure
ax.set_xlim(min(x), max(x))
ax.set_ylim(min(y), max(y))
ax.set_zlim(0, hauteur)
# Afficher l'image
plt.show()
else:
print("Aucun bâtiment trouvé à proximité de ces coordonnées.")
Mais ce n’est pas tout! Pourquoi ne pas extraire les 10 bâtiments les plus proches de notre bâtiment sélectionné ? Ça vous plaît, n’est-ce pas ? Voici comment on va procéder…
[ ]
# Trouver les 10 bâtiments les plus proches du point spécifié
dix_plus_proches = GDF.distance(point_batiment).nsmallest(11).index
batiments_proches = GDF.loc[dix_plus_proches]
batiments_proches_utm = batiments_proches.to_crs(epsg=26191)
# Vérifier si les bâtiments sont trouvés
if not batiments_proches_utm.empty:
# Créer une figure pour la visualisation
fig, ax = plt.subplots()
# Tracer les contours des bâtiments proches
batiments_proches_utm.plot(ax=ax, color='blue', alpha=0.5)
# Mettre en évidence le bâtiment spécifié
batiment2.plot(ax=ax, color='red')
# Ajouter un point pour les coordonnées du bâtiment spécifié
#ax.scatter(longitude_batiment, latitude_batiment, c='red', marker='x')
# Afficher l'image
plt.show()
else:
print("Aucun bâtiment trouvé à proximité de ces coordonnées.")
Et comme si ça ne suffisait pas, nous allons rajouter les données du réseau routier autour du bâtiment seul, puis autour du bâtiment et des 10 bâtiments les plus proches…
[ ]
# Assurez-vous que 'gdf_edges' utilise la même projection que 'GDF2'
gdf_edges = gdf_edges.to_crs(epsg=26191)
# Utiliser la géométrie du bâtiment sélectionné pour définir la zone d'intérêt
zone_interet = batiment2.unary_union.buffer(50) # Ajustez la valeur du buffer selon vos besoins
# Filtrez les arêtes du réseau routier qui intersectent la zone d'intérêt
edges_within_zone = gdf_edges[gdf_edges.intersects(zone_interet)]
if not edges_within_zone.empty:
fig, ax = plt.subplots()
batiment2.plot(ax=ax, color='red')
# Tracer les arêtes du réseau routier sur la figure
edges_within_zone.plot(ax=ax, color='grey', alpha=0.5, linewidth=1)
# Afficher l'image avec les bâtiments et le réseau routier pertinent
plt.show()
else:
print("Aucun réseau routier trouvé dans la zone d'intérêt.")
[ ]
# Assurez-vous que 'gdf_edges' utilise la même projection que 'GDF2'
gdf_edges = gdf_edges.to_crs(epsg=26191)
# Utiliser la géométrie du bâtiment sélectionné pour définir la zone d'intérêt
zone_interet = batiments_proches_utm.unary_union.buffer(50) # Ajustez la valeur du buffer selon vos besoins
# Filtrez les arêtes du réseau routier qui intersectent la zone d'intérêt
edges_within_zone = gdf_edges[gdf_edges.intersects(zone_interet)]
if not edges_within_zone.empty:
fig, ax = plt.subplots()
batiment2.plot(ax=ax, color='red')
batiments_proches_utm.plot(ax=ax, color='blue', alpha=0.5)
# Tracer les arêtes du réseau routier sur la figure
edges_within_zone.plot(ax=ax, color='grey', alpha=0.5, linewidth=1)
# Afficher l'image avec les bâtiments et le réseau routier pertinent
plt.show()
else:
print("Aucun réseau routier trouvé dans la zone d'intérêt.")
Vous voulez convertir ce dessin en un fichier DXF compatible avec AutoCAD ? Pas de problème…
[ ]
# Créer un fichier DXF
doc = ezdxf.new(dxfversion='R2010')
msp = doc.modelspace()
# Ajouter les bâtiments au DXF
for _, building in batiments_proches_utm.iterrows():
# Supposons que 'geometry' est une colonne dans votre géo-dataframe qui contient les polygones des bâtiments
if building['geometry'].geom_type == 'Polygon':
points = list(building['geometry'].exterior.coords)
msp.add_lwpolyline(points, close=True, dxfattribs={'layer': 'Bâtiments'})
elif building['geometry'].geom_type == 'MultiPolygon':
for poly in building['geometry']:
points = list(poly.exterior.coords)
msp.add_lwpolyline(points, close=True, dxfattribs={'layer': 'Bâtiments'})
# Ajouter les rues au DXF
for _, street in edges_within_zone.iterrows():
# Supposons que 'geometry' est une colonne dans votre géo-dataframe qui contient les lignes des rues
if street['geometry'].geom_type == 'LineString':
points = list(street['geometry'].coords)
msp.add_lwpolyline(points, dxfattribs={'layer': 'Rues'})
elif street['geometry'].geom_type == 'MultiLineString':
for line in street['geometry']:
points = list(line.coords)
msp.add_lwpolyline(points, dxfattribs={'layer': 'Rues'})
# Enregistrer le fichier DXF
doc.saveas('plan_bati.dxf')
print('Le fichier DXF a été généré avec succès.')
Le fichier DXF a été généré avec succès.
Et enfin, juste pour le plaisir, revenons à notre bâtiment et visualisons-le en 3D, dans une forme tridimensionnelle dynamique cette fois-ci, sous la projection EPSG26191, puis nous allons créer un fichier PLY à partir des données pour une utilisation ultérieure dans un logiciel de modélisation 3D.
[ ]
# Vérifier si le bâtiment est trouvé
if not batiment2.empty:
# Obtenir les coordonnées x, y, z du polygone
coords = batiment2.geometry.iloc[0].exterior.coords
x, y = zip(*coords)
z = np.zeros(len(x)) # Bas du bâtiment
# Spécifier une hauteur pour le bâtiment
hauteur = 50 # Hauteur en mètres, remplacez par la hauteur réelle si disponible
# Créer les points supérieurs du bâtiment en ajoutant la hauteur
z_top = np.full(len(x), hauteur)
# Définir les listes pour les indices des sommets de chaque face
I, J, K = [], [], []
# Générer les indices des sommets pour les faces verticales
for i in range(len(x)-1):
I.extend([i, i+1, i+1, i])
J.extend([i+1, i, i+len(x), i+len(x)])
K.extend([i+len(x), i+len(x)+1, i+1+len(x), i+1])
# Ajouter les indices pour le toit et le sol
# Indices pour le toit
I.extend([i for i in range(len(x))])
J.extend([i+len(x) for i in range(len(x))])
K.extend([i+1 if i < len(x)-1 else 0 for i in range(len(x))])
# Indices pour le sol
I.extend([i+len(x) for i in range(len(x))])
J.extend([i+1+len(x) if i < len(x)-1 else len(x) for i in range(len(x))])
K.extend([i+1 if i < len(x)-1 else 0 for i in range(len(x))])
# Créer la figure 3D interactive avec plotly
fig = go.Figure(data=[
go.Mesh3d(
x=np.concatenate([x, x]), # x-coordonnées
y=np.concatenate([y, y]), # y-coordonnées
z=np.concatenate([z, z_top]), # z-coordonnées
i=I, # Indices des sommets de chaque face (I)
j=J, # Indices des sommets de chaque face (J)
k=K, # Indices des sommets de chaque face (K)
facecolor=['rgba(0, 255, 255, 0.5)'], # Couleur des faces
opacity=0.5 # Transparence
)
])
# Mise à jour des paramètres de la figure
fig.update_layout(
scene=dict(
xaxis=dict(nticks=4, range=[min(x), max(x)]),
yaxis=dict(nticks=4, range=[min(y), max(y)]),
zaxis=dict(nticks=4, range=[0, hauteur])
),
width=700,
margin=dict(r=20, l=10, b=10, t=10)
)
# Afficher la figure
fig.show()
else:
print("Aucun bâtiment trouvé à proximité de ces coordonnées.")
def export_plotly_to_ply(fig, filename):
# Récupérer les données de la figure Plotly
mesh = fig.data[0]
x, y, z = mesh.x, mesh.y, mesh.z
i, j, k = mesh.i, mesh.j, mesh.k
# Créer l'en-tête du fichier PLY
ply_content = "plyn"
ply_content += "format ascii 1.0n"
ply_content += "element vertex {}n".format(len(x))
ply_content += "property float xn"
ply_content += "property float yn"
ply_content += "property float zn"
ply_content += "element face {}n".format(len(i))
ply_content += "property list uchar int vertex_indicesn"
ply_content += "end_headern"
# Ajouter les sommets au fichier PLY
for vx, vy, vz in zip(x, y, z):
ply_content += "{} {} {}n".format(vx, vy, vz)
# Ajouter les faces au fichier PLY
for vi, vj, vk in zip(i, j, k):
ply_content += "3 {} {} {}n".format(vi, vj, vk)
# Écrire le contenu dans un fichier PLY
with open(filename, "w") as file:
file.write(ply_content)
print(f"Les données ont été exportées dans '{filename}'.")
# Utiliser la fonction pour exporter la figure Plotly en fichier PLY
export_plotly_to_ply(fig, "batiment.ply")
Et qu’en est-il du traitement de données d’élévation ? Dans la cellule suivante, nous allons charger une image géoréférencée TIF, que nous possédons, représentant les données d’élévation d’une zone de la ville de Fès, y compris la zone d’intérêt. Nous allons ensuite extraire les données d’élévation depuis l’image pour créer une représentation 3D de la zone couverte dans l’image.
[ ]
!wget -O ASTGTM2_N34W006_dem.tif "https://ellipsoide.xyz/pypro/ASTGTM2_N34W006_dem.tif"
[ ]
# Charger l'image géoréférencée (remplacez 'votre_image.tif' par le chemin vers votre fichier GeoTIFF)
chemin_image = '/content/ASTGTM2_N34W006_dem.tif'
with rasterio.open(chemin_image) as src:
elevation_data = src.read(1) # Lire les données d'élévation (bande 1)
transform = src.transform # Obtenir la transformation géospatiale
# Créer une grille de coordonnées géographiques
xel, yel = np.meshgrid(np.arange(src.width), np.arange(src.height))
lonel, latel = rasterio.transform.xy(transform, yel, xel) # Convertir les coordonnées de pixel en coordonnées géographiques
# Afficher l'image
plt.figure(figsize=(8, 8))
plt.imshow(elevation_data, cmap='terrain')
plt.title('Image d'élévation')
plt.axis('off') # Masquer les axes
plt.show()
[ ]
# Créer une figure 3D
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# Afficher le terrain en 3D
surf = ax.plot_surface(lonel, latel, elevation_data, cmap='terrain', edgecolor='none')
# Ajouter des étiquettes
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_zlabel('Élévation')
# Ajouter un contrôle interactif pour zoomer et faire pivoter
#ax.view_init(elev=30, azim=-45) # Angle de vue initial
fig.colorbar(surf, shrink=0.5, aspect=10) # Ajouter une barre de couleur
# Afficher la figure
plt.title('Graphique 3D interactif de l'élévation du terrain')
plt.show()
Et pour situer notre bâtiment sur le terrain et extraire les données d’élévation associées aux coordonnées du bâtiment sélectionné précédemment, nous allons convertir les coordonnées géographiques du bâtiment en coordonnées de pixel puis nous allons recadrer notre zone d’intérêt pour situer le bâtiment sur le terrain et détecter la valeur d’élévation.
[ ]
# Convertir les coordonnées géographiques en coordonnées de pixel
pixel_x, pixel_y = rasterio.transform.rowcol(transform, longitude_batiment, latitude_batiment)
# Vérifier si le pixel est à l'intérieur de l'image
if 0 <= pixel_x < src.width and 0 <= pixel_y < src.height:
# Extraire l'élévation à ce pixel
elevation_at_point = elevation_data[pixel_y, pixel_x]
print(f"Élévation au point ({latitude_batiment:.6f}, {longitude_batiment:.6f}) : {elevation_at_point} mètres")
else:
print("Le point se trouve en dehors de l'image géoréférencée.")
# Définir la région d'intérêt (périmètre de 3000 mètres autour du point central)
min_x, max_x = pixel_x - 65, pixel_x + 65
min_y, max_y = pixel_y - 65, pixel_y + 65
# Créer une figure 3D pour la région d'intérêt
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# Afficher le terrain en 3D pour la région d'intérêt
surf_bat = ax.plot_surface(lonel[min_y:max_y], latel[min_y:max_y], elevation_data[min_y:max_y], cmap='terrain', edgecolor='none')
# Ajouter des étiquettes
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_zlabel('Élévation')
# Ajouter un contrôle interactif pour zoomer et faire pivoter
ax.view_init(elev=30, azim=-45) # Angle de vue initial
fig.colorbar(surf_bat, shrink=0.5, aspect=10) # Ajouter une barre de couleur
# Ajouter un point au graphique à la position du bâtiment
ax.scatter(longitude_batiment, latitude_batiment, elevation_at_point, color='red', s=100, label='Bâtiment', marker='x', zorder=9999999)
# Afficher la légende
ax.legend()
# Afficher la figure
plt.title(f"Élévation au point ({latitude_batiment:.6f}, {longitude_batiment:.6f}) : {elevation_at_point} mètres")
plt.show()
À présent, en utilisant les mêmes données d’élévation extraites depuis l’image, nous allons créer une figure pour représenter les courbes de niveau, et une seconde figure pour représenter le profil longitudinal.
[ ]
# Définissez le chemin vers l'image satellite de Fès (remplacez par le chemin réel)
image_path = '/content/ASTGTM2_N34W006_dem.tif'
# Ouvrez l'image avec Rasterio
with rasterio.open(image_path) as src:
# Lisez la bande d'élévation (par exemple, bande 1)
elevation = src.read(1).astype(float)
# Définir les niveaux de contour pour les courbes de niveau
contour_levels = np.arange(start=np.min(elevation), stop=np.max(elevation), step=50)
# Créer une figure pour les courbes de niveau
plt.figure(figsize=(8, 8))
plt.title("Courbes de niveau de l'image satellite de Fès")
plt.imshow(elevation, cmap='terrain') # Utilisez une colormap appropriée pour l'élévation
plt.contour(elevation, levels=contour_levels, colors='black', alpha=0.5)
plt.colorbar(label='Élévation (mètres)')
plt.axis('off')
plt.show()
[ ]
# Définissez le chemin vers l'image satellite de Fès (remplacez par le chemin réel)
image_path = '/content/ASTGTM2_N34W006_dem.tif'
# Ouvrez l'image avec Rasterio
with rasterio.open(image_path) as src:
# Lisez la b1nde d'élévation (par exemple, bande 1)
elevation = src.read(1).astype(float)
# Sélectionnez une ligne pour le profil longitudinal (par exemple, une colonne centrale)
row_index = elevation.shape[0] // 2
elevation_profile = elevation[row_index, :]
# Créez un graphique du profil longitudinal
plt.figure(figsize=(8, 4))
plt.plot(elevation_profile, label="Profil longitudinal")
plt.xlabel("Distance (pixels)")
plt.ylabel("Élévation (mètres)")
plt.title("Profil longitudinal de l'image satellite de Fès")
plt.grid()
plt.legend()
plt.show()
Dans la section suivante, nous allons charger, traiter, analyser et visualiser les données depuis deux images satellites de la ville de Fès . La première image contient les trois bandes RGB (Rouge, Vert, Bleu). La deuxième image contient la bande 3 (Rouge) ainsi que les bandes 6 et 7 situées dans la région du spectre infrarouge à ondes courtes (SWIR).
Donc nous commencerons d’abord par charger et afficher la première image satellite, puis nous utiliserons la bande rouge et la bande bleue pour calculer l’indice de végétation NDVI, au point d’intérêt (bâtiment sélectionné précédemment), sur un graphique géoréférencé.
[ ]
!wget -O satrgb.tif "https://ellipsoide.xyz/pypro/satrgb.tif"
arrow_upwardarrow_downwardlinksettingsdeletemore_vert
[ ]
# Définissez le chemin vers l'image satellite de Fès (remplacez par le chemin réel)
image_source = '/content/satrgb.tif'
# Ouvrez l'image avec Rasterio
with rasterio.open(image_source) as src:
# Lisez les bandes de l'image (par exemple, bandes 4, 3 et 2 pour une image RVB)
band1 = src.read(1)
band2 = src.read(2)
band3 = src.read(3)
# Empilez les bandes pour créer une image RVB
datasat = np.dstack((band1, band2, band3))
# Affichez l'image
plt.figure(figsize=(8, 8))
plt.imshow(datasat)
plt.title("Image satellite de Fès")
plt.axis('off') # Masquez les axes
plt.show()
# Lisez les bandes rouge (B3) et proche infrarouge (B1)
B3 = src.read(3).astype(float)
B2 = src.read(2).astype(float)
B1 = src.read(1).astype(float)
transdata = src.transform
# Calculez l'Indice de Végétation Normalisé (NDVI)
ndvi = (B3 - B1) / (B3 + B1)
# Convertir les coordonnées géographiques en coordonnées de pixel
bati_x, bati_y = rasterio.transform.rowcol(transdata, longitude_batiment, latitude_batiment)
# Vérifier si le pixel est à l'intérieur de l'image
# Extraire l'élévation à ce pixel
batisat_point = ndvi[bati_y, bati_x]
print(f"Indice de végétation au point ({latitude_batiment:.6f}, {longitude_batiment:.6f}) : {batisat_point} (NDVI)")
# Obtenez les coordonnées géographiques réelles
left, bottom, right, top = src.bounds
# Affichez l'image NDVI avec les coordonnées géographiques réelles
plt.figure(figsize=(8, 8))
plt.imshow(ndvi, cmap='RdYlGn', extent=(left, right, bottom, top), origin='upper')
plt.title("Indice de Végétation (NDVI)")
plt.colorbar(label="Valeur NDVI")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
# Tracez l'emplacement du bâtiment sur l'image NDVI
plt.scatter(longitude_batiment, latitude_batiment, c='black', s=100, marker='X', label='Emplacement du bâtiment')
plt.legend()
plt.show()
Puis, cette fois-ci, nous allons charger afficher la deuxième image satellite et nous utiliserons les bandes infrarouges à ondes courtes 6 et 7 pour calculer l’indice des zones bâties NDBI au point d’intérêt de la même manière.
[ ]
!wget -O sat367.tif "https://ellipsoide.xyz/pypro/sat367.tif"
[ ]
# Définissez le chemin vers l'image satellite de Fès (remplacez par le chemin réel)
image_path = '/content/sat367.tif'
# Ouvrez l'image avec Rasterio
with rasterio.open(image_path) as src:
# Lisez les bandes de l'image
band7 = src.read(1)
band6 = src.read(2)
band3r = src.read(3)
# Empilez les bandes pour créer une image RVB
datasatb = np.dstack((band7, band6, band3r))
B7 = src.read(1).astype(float)
B6 = src.read(2).astype(float)
ndbi = (B6 - B7) / (B6 + B7)
transbativ = src.transform
# Convertir les coordonnées géographiques en coordonnées de pixel
batib_x, batib_y = rasterio.transform.rowcol(transbativ, longitude_batiment, latitude_batiment)
# Vérifier si le pixel est à l'intérieur de l'image
# Extraire l'indice à ce pixel
batiib_point = ndbi[batib_y, batib_x]
print(f"Indice de zones bâties au point ({latitude_batiment:.6f}, {longitude_batiment:.6f}) : {batiib_point} (NDBI)")
# Obtenez les coordonnées géographiques
left, bottom, right, top = src.bounds
# Affichez l'image RVB
plt.figure(figsize=(8, 8))
plt.imshow(datasatb)
plt.title("Image satellite bandes 3,6 et 7 de Fès")
plt.axis('off') # Masquez les axes
plt.show()
# Affichez le NDBI (Normalized Difference Built-up Index)
plt.figure(figsize=(8, 8))
plt.imshow(ndbi, cmap='RdYlGn', extent=(left, right, bottom, top), origin='upper')
plt.colorbar()
plt.title("Indice de zones bâties NDBI")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
# Tracez l'emplacement du bâtiment sur l'image NDBI
plt.scatter(longitude_batiment, latitude_batiment, c='black', s=100, marker='X', label='Emplacement du bâtiment')
plt.legend()
plt.show()
Et enfin, une petite subtilité dans le cadre de traitement des images satellites, nous allons utiliser la bande 2, contenue dans la première image, et la bande 6, contenue dans la deuxième image, pour calculer la densité d’eau NDWI au point d’intérêt.
[ ]
# Calculez la densité d'eau (NDWI)
ndwi = (B2 - B6) / (B2 + B6)
# Convertir les coordonnées géographiques en coordonnées de pixel
bati_x, bati_y = rasterio.transform.rowcol(transdata, longitude_batiment, latitude_batiment)
# Vérifier si le pixel est à l'intérieur de l'image
# Extraire l'élévation à ce pixel
batiwat_point = ndwi[bati_y, bati_x]
print(f"La densité d'eau au point ({latitude_batiment:.6f}, {longitude_batiment:.6f}) : {batiwat_point} (NDWI)")
# Obtenez les coordonnées géographiques réelles
left, bottom, right, top = src.bounds
# Affichez l'image NDWI avec les coordonnées géographiques réelles
plt.figure(figsize=(8, 8))
plt.imshow(ndwi, cmap='RdYlGn', extent=(left, right, bottom, top), origin='upper')
plt.title("Indice de densité d'eau (NDWI)")
plt.colorbar(label="Valeur NDWI")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
# Tracez l'emplacement du bâtiment sur l'image NDWI
plt.scatter(longitude_batiment, latitude_batiment, c='black', s=100, marker='X', label='Emplacement du bâtiment')
plt.legend()
plt.figure(figsize=(12, 6))
plt.plot(B2[bati_y, :], label='Bande bleue')
plt.plot(B6[bati_y, :], label='Bande infrarouge')
plt.xlabel("Position du pixel")
plt.ylabel("Réflectance")
plt.title("Profil Longitudinal")
plt.legend()
plt.show()
Et voici un graphique 3D superposant toutes les bandes utilisées par niveau de réflectance.
[ ]
with rasterio.open(image_path) as src:
dbands = src.read()
# Obtenir les coordonnées géographiques réelles
transform = src.transform
x_coords, y_coords = np.meshgrid(np.arange(dbands.shape[2]), np.arange(dbands.shape[1]))
lon, lat = transform * (x_coords, y_coords)
max_reflectance = np.maximum.reduce([band1, band2, band3, band6, band7])
colors = np.where(max_reflectance == band1, 'green', # Vert
np.where(max_reflectance == band2, 'blue', # Bleu
np.where(max_reflectance == band3, 'magenta', # Rouge
np.where(max_reflectance == band6, 'red', # Infrarouge 6
np.where(max_reflectance == band7, 'purple', 'transparent' ))))) # Infrarouge 7
colors = colors.reshape(-1)
# Créer le graphique tridimensionnel
fig = plt.figure(figsize=(12, 8))
ax_combined = fig.add_subplot(111, projection='3d')
# Superposer les coordonnées géographiques pour chaque bande
ax_combined.scatter(lon, lat, max_reflectance, c=colors, label='Dominance', s=1)
# Personnalisation du graphique
ax_combined.set_title('Superposition des bandes avec réflectance maximale')
ax_combined.set_xlabel('Longitude')
ax_combined.set_ylabel('Latitude')
ax_combined.set_zlabel('Réflectance maximale')
ax_combined.legend()
plt.tight_layout()
plt.show()
Voilà ! On va s’arrêter là, pour l’instant. J’espère que vous avez apprécié cette démonstration. Créer de la valeur nécessite du temps et un effort considérable, n’hésitez pas à me rejoindre sur Linkedin ou sur notre site web, je me ferai un plaisir de recevoir vos remarques et vos impressions.
Et ce n’est pas fini ! L’aventure pourrait reprendre de plus belle à n’importe quel moment. Tout dépendra de la suite des événements 🙂
Pour un accès direct au code et pour que vous puissiez le modifier et l’utiliser à votre guise, vous pouvez consulter Ce lien
Mr. Ali OUFRID
Ingénieur Topographe et Géomètre Expert.
Une référence dans le domaine de la topographie et de la cartographie au Maroc et aux nations unies.
Ellipsoide
Contactez notre Bureau d'Etudes et Travaux Topographiques et Cartographiques