2011-02-18 Petits exercices avec Gnuplot

http://www.k1ka.be/pics/stations.jpeg

Gnuplot est connu comme un outil puissant quoique rébarbatif, essentiellement pour sa capacité à tracer des fonctions mathématiques ou des graphiques de type conventionnel, comme les histogrammes, les nuages de points, etc.

On peut cependant l'utiliser pour des travaux plus diversifiés. Des données associées à des coordonnées géographiques peuvent par exemple tracer d'assez jolies cartes, comme celle figurant en haut de page. Elles n'ont pas la prétention de rivaliser avec celles produites par des outils dédiés à ce domaine (en particulier parce que ces cartes ne relèvent que d'une projection élémentaire dite, si je ne me trompe, «plate carrée»), mais elles permettent au moins d'appréhender des données de façon intuitive par rapport à leur contexte spatial.

Dans l'exemple qui suit, j'ai cherché à représenter les observations météorologiques connues sous le nom de METAR, qui proviennent pour la plupart d'aéroports et sont disponibles en permanence sur le site de la National Oceanic and Atmospheric Administration , plus précisément à cette adresse: http://weather.noaa.gov/weather/metar.shtml

Ces données sont dans un format ASCII à première vue assez obscur:

 2011/02/17 11:50
 EBBR 171150Z 07005KT 2100 BR FEW004 BKN006 04/03 Q1004 BECMG 3500 BR SCT006

et des librairies en divers langages (dont Perl) permettent de faciliter leur interprétation1.

Elles sont disponibles soit sous forme de données isolées propres à une seule station, soit sous forme de fichiers assez important rassemblant toutes les obsrvations d'une tranche horaire (cycles files)2.

Une première chose à considérer est le format sous lequel transmettre des données à gnuplot. Celui-ci attend le plus souvent des données sous forme de lignes dont les éléments sont séparés par des espaces. Par exemple

  # Latitude Longitude	Température
    -22.6    63.97      1
    -21.9    64.13      2
      4.47   51.2       5
      5.47   51.17      6
   (etc.)

Mais dès qu'on envisage d'effectuer des sélections parmi des volumes importants de données, il devient intéressant de les transmettre de manière dynamique plutôt qu'à partir d'un fichier statique. Gnuplot peut lire en entrée le résultat d'une commande externe, ce qui permet d'utiliser toutes les combinaisons des outils en ligne de commande chaînés les uns aux autres, tels que sed, cut, etc.

Mais toute cette tuyauterie me semble plutôt réservée à des utilisateurs du shell assez aguerris. Les exemples de la documentation 3 proposent assez souvent d'utiliser awk, langage que je n'ai jamais eu le courage d'aborder. J'ai trouvé beaucoup plus commode d'utiliser SQL, et en particulier sqlite qui permet de conserver ses données dans le répertoire de travail sans même avoir à maintenir un serveur de base de données.

Les données METAR ne contiennent pas de coordonnées géographiques, mais les aéroports auxquels elles sont associées sont identifiés par un code unique ICAO. Celui-ci pourra être utilisé comme index dans les tables. J'ai retrouvé un fichier rassemblant toutes ces coordonnées et importé celles-ci ainsi que quelques données éventuellement exploitables (la nationalité, l'altitude, etc.) dans une table sqlite3. Les scripts et les détails de ces opérations sont donnés de manière extensive dans une autre page pour ne pas alourdir celle-ci.

À la suite de quoi, il était déjà possible de sélectionner les stations en fonction de leur nationalité, soit de manière interactive:

  gv@spirou:~$ sqlite3 airports.db
  SQLite version 3.7.3
  Enter ".help" for instructions
  Enter SQL statements terminated with a ";"
  sqlite> .schema airports
  CREATE TABLE airports (
    ICAO varchar(4) primary key,
    COUNTRY varchar(2),
    LAT decimal(5,2),
    LON decimal(5,2),
    ALTITUDE mediumint,
    NAME varchar(16)
  );
  sqlite> .separator '  '
  sqlite> SELECT ICAO,NAME FROM airports WHERE COUNTRY LIKE 'BE' ORDER BY NAME;
  EBAW  ANTWERP/DEURNE
  EBBE  BEAUVECHAIN (BAF
  EBBX  BERTRIX (BEL-AFB
  EBLG  BIERSET/LIEGE (C
  ...
  EBSZ  SEMMERZAKE (MIL)
  EBSP  SPA/LA SAUVENIER
  EBST  ST. TRUIDEN (BAF
  EBWE  WEELDE (MIL)
  sqlite> .quit

soit depuis le shell. Combien de stations ICAO sont situées en France ?

  gv@spirou:~$ sqlite3 airports.db \
   "SELECT COUNT(ICAO) FROM airports WHERE COUNTRY LIKE 'FR'"
  149

Y a-t-il des aéroports sous le niveau de la mer ?

  gv@spirou:~$ sqlite3 -separator ' ' airports.db
  "SELECT ICAO,COUNTRY,NAME,ALTITUDE FROM airports WHERE ALTITUDE < 0"
  KNJK US EL CENTRO NAF    -13
  KIPL US IMPERIAL         -17
  KTRM US THERMAL/PALM SPG -35
  SKSP CO SAN ANDRES ISLAN -2
  EHLE NL LELYSTAD         -6
  EHRD NL ROTTERDAM AIRPOR -5
  OIGG IR RASHT            -13
  OINN IR NOSHAHR          -19
  OINR IR RAMSAR           -20
  UATG KZ ATYRAU           -22
  UBBB AZ BAKU/BINE ARPT   -6
  URWA RU ASTRAKHAN        -20

Apparemment, oui :) !

Voilà donc finie la première étape qui consistait à pouvoir localiser les données.

Pour décoder les données METAR proprement dite, j'ai utilisé le module Geo-METAR-1.15 de Koos van den Hout: http://search.cpan.org/~koos/Geo-METAR-1.15/METAR.pm.

Ces données seront insérées dans une seconde table, qui aura en commun avec la première le champs ICAO identifiant la station. Le problème est qu'elles sont parfois incomplètes ou mal formées. Or Gnuplot interrompra le tracé dès qu'il rencontrera une donnée erronée. Il faut donc veiller à vérifier la validité des données extraites avant de les insérer, et à ce la table ne contienne que des valeurs 'NULL' dans le cas contraire, ce qui permet de faire appel à cette condition dans les requêtes.

Au final, j'ai donc pu utiliser la requête suivante pour obtenir les coordonnées et la température des stations comprise dans une zone correspondant à l'Europe sans la Scandinavie (celle qui est représentée dans la petite carte en haut de page).

 $ sqlite3 -separator ' ' airports.db 'SELECT LON,LAT,TEMP FROM airports \
   LEFT JOIN metar ON metar.ICAO=airports.ICAO \
   WHERE lat < 65 and lat > 33 and lon > -25 and lon < 40 \
   AND TEMP IS NOT NULL '
 
   -22.6 63.97 1
   -21.9 64.13 2
   4.47 51.2 5
   5.47 51.17 6
   4.5 50.9 6
   4.45 50.47 6
   [...]

Pour en faire une carte compréhensible, il faudra y superposer le tracé des lignes côtières. Du moment qu'on dispose de ce type de données sous forme de couples de coordonnées, on peut les tracer avec gnuplot. Tant qu'on se contente d'une définition assez faible, on peut trouver assez facilement des données de ce type librement utilisables. J'en ai répertorié sur cette page, et pour cet exemple j'ai retenu celles du NOAA: http://www.ngdc.noaa.gov/mgg/coast/ qu'on peut sélectionner en précisant les latitudes et longitudes de la zone désirée.

Une possibilité aurait été de tracer un graphe en 2 dimensions en utilisant la valeur de la troisième colonne comme marqueur (c'est l'option with labels).

  plot "< sqlite3 -separator ' ' airports.db \
  'SELECT LON,LAT,TEMP FROM airports \
  LEFT JOIN metar ON metar.ICAO=airports.ICAO \
  WHERE lat < 65 and lat > 33 and lon > -25 and lon < 40 AND TEMP IS NOT NULL ' "\
  with labels  notitle,\
  'coastlines.dat' with lines lc rgb "grey" notitle

http://www.k1ka.be/pics/raw0.gif

Mais comme on le voit, il y a un problème de collision entre les températures...

Comme les données retournées sont des triplets, cela permet de les positionner dans un espace en 3 dimensions et d'utiliser la commande splot.

Pour arriver au résultat souhaité, nous utiliserons tout d'abord une fonctionnalité de Gnuplot qui permet de reconstituer une surface sous forme de filet à partir de points répartis irrégulièrement. Il s'agit de la fonction dgrid3d qui demande comme arguments 2 premiers entiers qui détermineront le nombre de mailles désiré selon les 2 axes horizontaux (et donc la finesse du résultat), et un troisième dont la définition est assez compliquée à résumer (je renvoie à la documentation de gnuplot). La hauteur de chaque point résultant est influencée par l'ensemble des points de donnée de manière inversement proportionnelle à leur distance, selon un calcul où intervient ce dernier facteur. En pratique, la surface s'adoucit quand celui-ci est plus grand.

C'est une fonction assez gourmande en ressources, et il vaut mieux respecter la recommandation de la doc consistant à utiliser des puissances de 2, ce qui accélère le rendu.

http://www.k1ka.be/pics/raw1.gif

Il est possible de représenter une telle surface en lui attribuant une couleur en fonction de la hauteur, c'est l'option set pm3d4, qui améliore beaucoup sa lisibilité; voici le résultat avec la palette par défaut et l'index des couleurs sur la droite.

 set xrange [-12:36]
 set yrange [33.75:62.5]
 set cbrange [-25:25]
 set dgrid3d 75,125,4
 unset surface
 set pm3d at s
 splot "< sqlite3 -separator ' ' airports.db 'SELECT LON,LAT,TEMP FROM \
 airports LEFT JOIN metar ON metar.ICAO=airports.ICAO \ 
 WHERE lat < 65 and lat > 33 and lon > -25 and lon < 40 \
 AND TEMP IS NOT NULL ' " \
 with pm3d at s  notitle

http://www.k1ka.be/pics/raw2.gif

C'est plus joli mais pas très explicite pour autant. Il faudra adopter un point de vue depuis le haut pour obtenir une carte. Il y a une instruction spécifique pour cela, c'est set view map.

http://www.k1ka.be/pics/raw3.gif

En ajoutant les lignes côtières, on arrive peu à peu au résultat souhaité:

http://www.k1ka.be/pics/temperatures.gif

Il est également possible de spécifier une palette de couleur différente. Je trouve que celle qui suit, en réservant les couleurs les plus froides à l'intervalle sous 0°, et les couleurs les plus claires aux extrémités de l'échelle, est assez explicite (ça se discute, mais je la préfère à la majorité de ce que j'ai trouvé sur le web5.

Enfin, j'y ai ajouté quelques lignes isothermes, reliant sur la surface obtenue tous les points de même hauteur, à la manière de courbes de niveaux.

http://www.k1ka.be/pics/temperatures2.gif

Le fichier de commande rassemblant toutes ces opérations est repris in extenso sur la page précitée, ainsi qu'un script Perl permettant de créer cette carte à la volée avec les dernières données disponibles.

Il s'y trouve aussi une carte animée pour un cycle de 24 heures, qui permet de réaliser que les températures nocturnes les plus froides n'ont pas lieu au milieu de la nuit mais plutôt avant le lever du soleil.

Mais en fin de compte, le résultat est-il fiable ? Le problème est que là où les points de coordonnées sont espacés, le graphe applique une «altitude» moyenne basée sur l'ensemble des données. Pour les zones couvertes par la mer, le résultat est vraissemblablement faussé à la fois par l'absence de relevés et parce que les conditions météorologiques à la surface de la mer sont très différentes de celles qu'on observerait dans des conditions égales au-dessus du sol. Il faut plutôt prendre ceci comme un exercice pour montrer les possibilités de Gnuplot; pour connaître le temps, rien de tel que de regarder par la fenêtre.

Je me demande par ailleurs si ce sont ces données METAR qui sont utilisées par tous les sites, commerciaux ou non, qui proposent des applets ou des bouts de code javascript permettant d'afficher le temps, sur le bureau ou une page web.

Je n'ai pas la réponse...

Footnotes:

1. Il existe aussi des décodeurs en ligne comme celui-ci: http://aviationweather.gov/adds/metars/index.php
3. Voir la section consacrée aux noms de fichiers particuliers http://www.gnuplot.info/docs_4.2/gnuplot.html#x1-13200034.1.7
5. comme ici par exemple)