Planète Rouge

Il y a un an je suis tombé sur une carte interactive de Mars (dont je ne retrouve pas l´adresse, forcément). Le résultat était très beau, mais pauvre en informations et non éditable. Comme l´article contenait quelques informations sur la manière dont le projet avait été réalisé, j´ai finalement décidé de faire ma propre carte - une bonne occasion de se frotter à des outils que je ne connaissais pas.

Générer les images

Les données

Pour construire ma carte, il me fallait avant toute chose des données. Par chance, les gens de la NASA sont sympa - les données du projet MOLA (Mars Orbiter Laser Altimeter) sont disponibles en ligne, pour tous ceux qui n´ont pas les moyens d´envoyer une sonde sur Mars. Je me suis donc retrouvé avec un gros fichier de données altimétriques allant de 88°N à 88°S, avec une résolution de 128 pixels par degré.

GDAL : le couteau suisse des GIS

A partir de ces données binaires, il fallait générer des images. J´ai donc appris les rudiments de l´utilisation de GDAL, une boite à outils permettant de manipuler les fichiers GIS (Geographic Information System).

Projection de Mercator

Avant de générer les images proprement dites, j´ai commencé par passer les données dans une projection adaptée à mon objectif : la projection de mercator.

gdalwarp --config GDAL_CACHEMAX 500 -wm 500 fichier_source.adf \
    -s_srs "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396190 +b=3396190 +units=m +no_defs" \
    -t_srs EPSG:3785 -r bilinear -te -20037508.34 -20037508.34 20037508.34 20037508.34 -ts ${RES} ${RES} \
    -co TILED=YES -co COMPRESS=lzw source_mercator.tif

Pour tout vous avouer, j´ai trouvé cette formule je ne sais plus où et je me suis abstenu de trop la modifier. ma seule contribution ici a été l´ajout du paramètre -ts ${RES} ${RES}, qui permet de redimensionner le résultat a une résolution donnée. Comme on le verra plus tard, il est préférable de travailler avec des images dont la résolution est une puissance de 256. Pour la version finale, j´ai opté pour RES=16384.

Cette étape passée, je pouvais me pencher sur la génération des images qui allaient peupler ma carte.

Rendu du relief

Pour commencer, j´ai défini une palette de couleurs pour représenter l´altitude. Mes critères de choix ont été surtout d´ordre esthétique. Un véritable cartographe hurlerait certainement d´horreur (d´autant plus qu´il n´y a même pas de légende sur ma carte !)

# altitude  rouge   vert    bleu
21249       250     227     158
0           250     144     66
-8208       194     77      20

J´ai ensuite appliqué ma palette de couleurs sur ma source mercator :

gdaldem color-relief -co COMPRESS=lzw source_mercator.tif ../palette_relief.txt image_relief.tif

Il faut savoir que l´image ainsi obtenue est assez «plate» quelque soit la palette choisie. pour faire ressortir les reliefs, il faut ombrer les escarpements. J´ai donc générer une seconde image :

gdaldem hillshade -combined -co COMPRESS=lzw source_mercator.tif image_hillshade.tif

Fusion d´images

J´ai longtemps galéré avant de trouver l´outil adapté pour fusionner les deux images. Impossible d´utiliser des outils comme Gimp ou ImageMagick sur des fichiers aussi lourds. J´ai finalement trouvé un tutoriel expliquant comment faire la fusion avec Mapnik.

Il me fallait donc un fichier de config …

<Map srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs">

    <Style name="color relief style">
        <Rule>
            <RasterSymbolizer mode="normal"/>
        </Rule>
    </Style>
    <Style name="hillshade style">
        <Rule>
            <RasterSymbolizer opacity="0.8" comp-op="multiply" />
        </Rule>
    </Style>

    <Layer name="color relief">
        <StyleName>color relief style</StyleName>
        <Datasource>
            <Parameter name="type">gdal</Parameter>
            <Parameter name="file">image_relief.tif</Parameter>
        </Datasource>
    </Layer>
    <Layer name="hillshade">
        <StyleName>hillshade style</StyleName>
        <Datasource>
            <Parameter name="type">gdal</Parameter>
            <Parameter name="file">image_hillshade.tif</Parameter>
        </Datasource>
    </Layer> 
</Map>

… et un script.

#!/usr/bin/python2
import mapnik
from sys import argv

res = int(argv[1])
print "resolution", res

map = mapnik.Map(res,res)
mapnik.load_map(map, 'mapnik_param.xml')
map.zoom_all()
mapnik.render_to_file(map, 'image_finale.tif')

Et me voilà enfin avec mon image finale, un TIFF de dimension 16384x16384, pesant 82Mo compressé.

Génération du tileset

Avec une carte aussi lourde (plus de 750Mo décompressée), il n´est bien entendu pas envisageable de tout charger d´un coup. C´est pourquoi les cartes interactives utilisent des tilesets. Un tileset est un ensemble d´images de dimension 256x256 dont l´assemblage forme l´image complête. Ainsi, la carte n´a besoin de charger que les tiles correspondant à la zone affichée. Pour gérer le zoom, il suffit de diviser la résolution de l´image par deux et de produire un nouveau tileset - et ce, autant de fois que nécessaire. Les tiles affichés seront alors choisi en fonction de la zone affichée et du zoom. C´est pour cette raison qu´il est préférable de partir d´une carte dont le dimensions sont une puissance de 256.

GDAL fournit un outil permettant de produire un tileset à partir d´une image.

gdal2tiles.py -p raster -s '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs' -z 0-6 -w none image_finale.tif surface

J´ai obtenu une arborescence contenant le tileset, qu´il ne restait plus qu´à afficher correctement sur une page web.

Et un deuxième tileset pour la route

En plus des images de terrain en fausses couleurs, j´ai aussi généré un deuxième tileset représentant les courbes de niveau. Là encore, GDAL s´occupe de tout.

gdal_contour -i 500 source_mercator.tif contour.shp

Cette commande génère un fichier de contours avec une ligne de niveau tous les 500 mètres. J´ai essayé avec un espacement de 100m et 200m, mais le résultat était bien trop chargé.

Après quoi, j´ai intégré les contours dans la carte d´ombrages pour accentuer les dénivelés.

Interface web

Il existe dans le monde merveilleux du logiciel libre un outil qui permet d´intégrer une carte interactive dans une page web en quelques lignes de javascript. Cet outil s´appelle Leaflet. Les tutoriels du site officiels sont très bien faits, je ne vais pas les répéter ici.

Tout l´intérêt du projet réside dans la possibilité d´afficher n´importe quel tileset, pour peu qu´il respecte une certaine nomenclature, et d´ajouter des éléments interactifs (marqueurs, zones, bulles de texte …)

Données géographiques

La carte étant construite, je me suis mis en quête d´informations à y intégrer - les noms des cratères, des plaines, des valées … Autant de choses que je n´avais pas l´intention de rentrer à la main.

J´ai finalement trouvé une base de données sur un serveur du MIT contenant les noms et coordonnées géographiques d´un grand nombres de particularités du terrain martien. Pour intégrer ces données dans ma carte, j´ai écrit un script python convertissant les données du format CSV au format GeoJSON utilisé par leaflet.

# !/usr/bin/python
# parse mars_db and produce GeoJson data
# http://stuff.mit.edu/afs/athena/project/xephem/lib/auxil/mars_db

import csv
import json

reader = csv.reader(open("mars_db.csv"),delimiter="|",skipinitialspace=True)

features = {}
        
for row in reader:
    name = row[0].strip()
    lat = float(row[1])
    lon = 360-float(row[2])
    dia = float(row[3])
    info = row[4].split(" ",1)[0].strip(",").lower()

    if info not in features:
        features[info] = []

    features[info].append({
        "type": "Feature",
        "properties": {
            "name": name,
            "type": info,
            "diameter": dia},
        "geometry": {
            "type": "Point",
            "coordinates": [lon,lat]}
        })


f = open("geo.json","w")

f.write("""var geo_data = {\n""")
for name in sorted(features.keys()):
    f.write(""""%s": { "type": "FeatureCollection", "features": %s },\n""" %(name,json.dumps(features[name],sort_keys=True)))
f.write("""}""");

f.close()

Affichage des données géographiques

Le nombre de marqueurs ainsi générés étant assez grand (plus de 500), j´ai décidé de faire un layer séparé pour chaque type de marqueur, avec la possibilité d´afficher ou non chaque layer.

![](/media/data.jpg

Ce qu´il reste à faire

Telle quelle, la carte est fonctionnelle, mais certains sont encore à ajouter :

  • remplacer les icones des marqueurs par quelque chose de plus esthétique ;
  • intégrer une fonction de recherche par texte ;
  • faire apparaitre automatiquement les bulles d´info au survol, en tenant compte du du champ diameter des marqueurs.