Affichage 3D des données OpenStreetMap sur une surface SRTM

Quoi ?

Afficher des chemins de randonnée issus d’OpenStreetMap en 3D pour avoir une meilleure idée du trajet dans la montagne.

Comment ?

Python et d’excellentes librairies (Mayavi mlab, gdal, Numpy, et Osm Python API) pour apprendre à manipuler ces données. C’est encore un travail en cours. J’ai utilisé pas mal de pages web/tutoriaux/documentations pour faire ce petit projet (voir liens).

J’ai dû ajouter le support des proxys à l’API Python d’OSM pour pouvoir l’utiliser.

J’ai choisi un chemin de randonnée près de Grenoble dans OpenStreetMap, cliqué sur « editer » pour obtenir l’ID du chemin. Le fichier SRTM correspondant à Grenoble (38_03) à été téléchargé ici pour obtenir les données d’élévation du terrain.

Prochaines étapes :

• télécharger et afficher plusieurs chemins, voir d’autres données (avec OsmApi, ça devrait être facile)
• obtenir de meilleures données d’élévation (difficile, car les données de l’IGN  à 25m de résolution ne sont pas libres)
• obtenir les traces GPS d’OpenStreetMap pour obtenir des informations d’élévation plus précises (nécéssite une mise à jour de l’API OSM Python)
• Afficher des textures sur le terrain, images satellites par exemple
• Afficher comme texture des rendus OSM existants (incluant forêt/roche/prairies…)

Voici le code (les commentaires devraient être assez clairs)

# -*- coding: utf-8 -*-
# ipython -wthread

from enthought.mayavi import mlab
import numpy
from osgeo import gdal, osr
import OsmApi
import math

earthRadius = 6371000 # rayon de la terre en mètres (oui c'est une approximation) https://en.wikipedia.org/wiki/Earth_radius

# lecture des données d'élévation SRTM
# Lecture du GeoTiff and récupération des coordonnées
ds = gdal.Open('/home/.../OpenStreetMap/SRTM/srtm_38_03.tif')

width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
maxx = gt[0] + width*gt[1] + height*gt[2]
miny = gt[3] + width*gt[4] + height*gt[5] 
maxy = gt[3]

band = ds.GetRasterBand(1)

# Zone de Grenoble qu'on veut afficher = pixels 400->1400, 5300->6000
zoneXmin = 400
zoneXsize = 1000
zoneYmin = 5300
zoneYsize = 695
data = band.ReadAsArray(400, 5300, 1000, 695).astype(numpy.float)

def degrees2metersLongX(latitude, longitudeSpan):
  """ latitude (in degrees) is used to convert a longitude angle to a distance in meters """
  return 2.0*math.pi*earthRadius*math.cos(math.radians(latitude))*longitudeSpan/360.0

def degrees2metersLatY(latitudeSpan):
  """ Convert a latitude angle span to a distance in meters """
  return 2.0*math.pi*earthRadius*latitudeSpan/360.0

def degrees2meters(longitude, latitude):
  return (degrees2metersLongX(latitude, longitude), degrees2metersLatY(latitude))

# Calcul des coordonnées X,Y pour le tableau data (qui contient l'altitude Z en mètres)
linex = numpy.arange(minx+zoneXmin*gt[1], minx+(zoneXmin+zoneXsize)*gt[1], gt[1])[0:zoneXsize]
arrayx = numpy.tile(linex, (len(data), 1)).transpose()

lineydegrees = numpy.arange( maxy+zoneYmin*gt[5], maxy+(zoneYmin+zoneYsize)*gt[5], gt[5])[0:zoneYsize]
liney = numpy.array([degrees2metersLatY(j) for j in lineydegrees])
arrayy = numpy.tile(liney, (len(data[0]), 1))

# Recalcule les position en mètres dans arrayx
for x, y in numpy.ndindex(arrayx.shape):
  arrayx[x,y] = degrees2metersLongX(lineydegrees[y], arrayx[x,y])

zscale = 1

# Affichage de la surface 3D
mlab.mesh(arrayx, arrayy, data.transpose())

# OpenStreetMap
osm = OsmApi.OsmApi(proxy="******", proxyport = 3128)

# Chemin du Bec du Corbeau : http://www.openstreetmap.org/browse/way/47189761 --> ID 47189761
way = osm.WayGet(47189761)
nodesId = way[u'nd']
nodes = osm.NodesGet(nodesId)
allX = []
allY = []
allZ = []

for nid in nodes:
  n=nodes[nid]
  (x,y) = degrees2meters(n['lon'], n['lat'])
  allX.append(x)
  allY.append(y)
  zz = data.transpose()[int(round((n['lon'] - (minx+zoneXmin*gt[1])) / gt[1])), int(round((n['lat'] - (maxy+zoneYmin*gt[5])) / gt[5]))]
  allZ.append(zz+30.0) # On affiche le chemin 30m au-dessus de la surface pour qu'il soit visible

# Affichage des noeuds du chemin sous forme de sphères
mlab.points3d(allX,allY,allZ, color=(1,0,0), mode='sphere', scale_factor=30)
# L'affichage d'une ligne ne fonctionne pas car les noeuds ne sont pas dans l'ordre logique
# mlab.plot3d(allX,allY,allZ, color=(1,1,0), line_width=15.0) # representation='surface' 'wireframe' 'points'

# Affichage d'une sphère noire aux  coordonnées de Grenoble
mlab.points3d([5.717*degrees2meters], [45.183*degrees2meters],[200], color=(0,0,0), mode='sphere', scale_factor=50)

Résultat

• Le « Y » des trois vallées de Grenoble

• Une vue zoomée pour mieux voir le chemin