Une carte météo avec Gnuplot

Gnuplot

Cette page explique comment les observations disponibles dans les rapports METAR (METeorological Airport Report) peuvent être visualisées sur une carte créée avec gnuplot.

Ces observations sont téléchargeables sur le site ftp du NOAA, le National Weather Service: ftp://tgftp.nws.noaa.gov/data/observations/metar/

Ces données ne comprennent pas les coordonnées géographiques des stations. Quelques recherches m'ont cependant permis de trouver ces renseignements dans un fichier texte disponible là: http://www.rap.ucar.edu/weather/surface/stations.txt Malheureusement, les noms des aéroports sont tronqués12.

Un peu de travail est cependant nécessaire pour pouvoir les exploiter:

  $ sqlite3 test.db
  SQLite version 3.6.21
  Enter ".help" for instructions
  Enter SQL statements terminated with a ";"
  sqlite> 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> .import 'icao.txt' airports
  sqlite> SELECT * FROM airports LIMIT 1;
  PADK:US:51.88:-176.65:4:ADAK NAS

Dans ce qui suit, on prendra comme limites géographiques une zone située entre les latitudes 33.75N et 62.5N, et les longitudes -12 et 36, ce qui couvre une bonne partie de l'Europe. Les lignes côtières seront tracées grâce à un fichier de coordonnées obtenu sur ce site:

http://www.ngdc.noaa.gov/mgg/coast/

On peut visualiser ces stations via ce fichier de commandes gnuplot:

  set term gif size 800,600
  set output 'stations.gif'
  set xrange [-12:36]
  set yrange [33.75:62.5]
  set grid
  plot "< sqlite3 -separator ' ' airports.db 'SELECT LON,LAT FROM airports \   
  WHERE lat < 62.5 and lat > 33.75 and lon > -12 and lon < 36' " \
  with points notitle, \
  'coastlines.dat' with lines lc rgb "red" notitle

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

Le fichier de commande final

  set term gif size 800,600
  set output 'output.gif'
  set multiplot
  set origin 0.5,0.5
  set xrange [-12:36]
  set yrange [33.75:62.5]
  set cbrange [-25:25]
  set dgrid3d 75,125,4
  set view map
  set tmargin at screen 0.95
  set bmargin at screen 0.05
  set lmargin at screen 0.05
  set rmargin at screen 0.85 # laisser de la place pour la palette
  unset surface
  set pm3d at b
  set grid
  set palette defined (0 0.8 1 1, 40 0.2 0.3 0.8, 60 0.3 0.8 0.2, 90 1 0.6 0, 100 1 1 0.2) 
  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 ' " notitle
  
  #-----------------------------
  set format
  set notics
  unset colorbox
  
  # tracé des isolignes
  set origin 0.5,0.5
  set contour base
  unset surface
  set cntrparam bspline
  set cntrparam levels discrete -10,0,10
  set style increment userstyle
  set style line 2 lc rgb "#808080"
  set style line 3 lc rgb "#b0b0b0"
  set style line 4 lc rgb "#505050"
  set grid
  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 lines notitle
  
  
  # tracé des lignes côtières
  set origin 0.5,0.5
  plot 'coastlines.dat' with lines lc rgb "white" notitle
  unset multiplot

Les scripts

extract_stations.pl

Le script qui a permis d'insérer dans la base données les coordonnées des stations.

  #!/usr/bin/perl
  
  open STATIONS, "stations.txt";
  while (<STATIONS>) {
    if (/^.{2} (.{16}) (\w{4}) .{14}(.{6})  (.{7}) (.{4})   X.*(..)$/) {
      $nom  = $1;
      $ICAO = $2;
      $lat  = $3;
      $lon  = $4;
      $alt  = $5;
      $pays = $6;
  
      $lat =~ /(..) (..)(.)/;
      $signe = $3 eq S ? '-' : '';
      $lat = $signe . sprintf("%.2f",$1+$2/60);
      $lon =~ /(...) (..)(.)/;
      $signe = $3 eq W ? '-' : '';
      $lon = $signe . sprintf("%.2f",$1+$2/60);
  
      $alt =~ s/^\s*//;
  
      print "$ICAO:$pays:$lat:$lon:$alt:$nom\n";
    }
  }

last_temp.pl

Un script permettant en une fois de

  #!/usr/bin/perl
  
  use strict;
  use LWP;
  use Geo::METAR;
  use DBI;
  
  my @gmtime = gmtime(time);
  my ($hour,$day, $month, $year) = (gmtime)[2,3,4,5];
  $hour = sprintf("%02d", $hour); # heure UTC
  my $date = sprintf("%04d-%02d-%02d", $year+1900, $month+1, $day);
  my $title = "$date - $hour h UTC";
  print "$title\n";
  my $file = $hour . 'Z.TXT';
  my $url = 'ftp://tgftp.nws.noaa.gov/data/observations/metar/cycles/' 
            . $file;
  my $browser = LWP::UserAgent->new();
  $browser->agent("libwww-perl 5.836-1");
  my $res = $browser->request(HTTP::Request->new(GET => $url));
  if ($res->is_success) {
      print $res->status_line, "\n";
      my $content=$res->content;
      &insert_in_db($content);
      &call_gnuplot
  } else {
      print $res->status_line, "\n";
  }
  
  sub call_gnuplot {
    my $commands=<<EoC;
    set term gif size 800,600
    set output '${hour}.gif'
    set title '$title'
    set multiplot
    set origin 0.5,0.5
    set xrange [-12:36]
    set yrange [33.75:62.5]
    set cbrange [-25:25]
    set dgrid3d 75,125,4
    set view map
    set pm3d at b
    set grid
    set tmargin at screen 0.90
    set bmargin at screen 0.05
    set lmargin at screen 0.05
    set rmargin at screen 0.85 # laisser de la place pour la palette
    unset surface
    set palette defined (0 0.8 1 1, 40 0.2 0.3 0.8, 60 0.3 0.8 0.2, 90 1 0.6 0, 100 1 1 0.2) 
    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 b  notitle
    
    #-----------------------------
    set format
    set notics
    unset colorbox
    
    # tracé des isolignes
    set origin 0.5,0.5
    set contour
    set cntrparam bspline
    set cntrparam levels discrete -30,-20,-10,0,10
    set style increment userstyle
    set style line 2 lc rgb "#808080"
    set style line 3 lc rgb "#b0b0b0"
    set style line 4 lc rgb "#505050"
    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 lines notitle
    
    # tracé des lignes côtières
    set origin 0.5,0.5
    plot 'coastlines.dat' with lines lc rgb "white" notitle
    unset multiplot
  EoC
    # fin du fichier de commandes gnuplot
    
    
    print "Calling Gnuplot\n";
    open GNUPLOT,"| gnuplot '-' 2> /dev/null";
    print GNUPLOT $commands;
    close GNUPLOT;
  }
  
  sub insert_in_db {
    my @lines = split(/\n/,$_[0]);
    # my $count = @lines; print "$count\n";
    my $metars={}; # création d'une table de hachage avec les noms ICAO comme
                   # clés, pour éviter les doublons dans la base sqlite.
    my $stations = 0; # compteur pour les stations
    my $temps    = 0; # compteur pour les températures
    my $dbh = DBI->connect( "dbi:SQLite:airports.db" ) 
              || die "Cannot connect: $DBI::errstr";
    print "Connecting to airports.db\n";
    $dbh->do( "DELETE FROM `metar`"); # on vide la table
    
    foreach (@lines) {
      if (/^\w{4} /) { # les lignes de données débutent par un nom ICAO
                       # composé de 4 lettres
        my $m = new Geo::METAR;
        $m->metar($_);
        my $site = $m->SITE;
        unless (defined $metars->{$site}) {
          $metars->{$site} = 1;
          $stations++;
          my $values = "'$site'"; my $columns = 'ICAO';
          my $temp_c = $m->TEMP_C;
          if ($temp_c =~ /^-*[0-9.]+$/) {
            $values = $values . ',' . "'$temp_c'";
            $columns = $columns . ',' . 'TEMP';
            $temps++;
          } 
          my $wind_dir_deg = $m->WIND_DIR_DEG;
          if ($wind_dir_deg =~ /^\d+$/) {
            $values = $values . ',' . "'$wind_dir_deg'";
            $columns = $columns . ',' . 'WIND_DIR_DEG';
          }
          my $wind_ms      = $m->WIND_MS;
          if ($wind_ms =~ /^\d+$/) {
            $values = $values . ',' . "'$wind_ms'";
            $columns = $columns . ',' . 'WIND_MS';
          }
          my $alt = $m->ALT;
          if ($alt =~ /^\d+$/) {
            $values = $values . ',' . "'$alt'";
            $columns = $columns . ',' . 'ALT';
          }
          # print "INSERT INTO metar ($columns) VALUES ($values)\n";
          $dbh->do( "INSERT INTO metar ($columns) VALUES ($values)");
        }
      }
    }
  print "  $stations stations processed.\n";
  print "  $temps temp. collected.\n"
  }
  

Une carte animée sur 24h

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

Footnotes:

1. Noter que sur Openstreetmap, les contributeurs semblent avoir enregistré le code ICAO des aéroports. Il devrait donc être possible d'obtenir ces données au format xml, en partant d'une requête de ce type
  http://xapi.openstreetmap.org/api/0.6/node[icao=*][bbox=-180,-90,180,90]

Mais le problème est que le taux de charge du serveur ne permet manifestement pas d'assumer une requête sur l'ensemble du monde.

2. Une autre ressource, au format csv: http://www.ourairports.com/data/