dimanche 3 novembre 2013

Mise en oeuvre des modèles météo WRF-ARW et WRF-NMM - Partie 6 - Visualisation des données

Dans le précédent article, vous aviez enfin, après de nombreuses et complexes manipulations, réussi à faire fonctionner votre premier run WRF-ARF. Il est grand temps de mettre un peu d'image derrière ces austères lignes de commande et de visualiser le résultat de notre travail. Aujourd'hui, je vous propose un mini-tutorial sur le langage de commande NCL qui permet de créer des scripts de visualisation scientifique sous forme de cartes ou de diagrammes de différents types.

Pour bien suivre ce tutorial

 
NCL est avant tout un outil en ligne de commande. Pour être utilisé convenablement, il convient de le mettre dans votre PATH et de paramétrer la variable d'environnement NCARG_ROOT. Si ce n'est pas déjà fait, je vous invite à consulter la partie 3 de ce tutorial où est expliquée la marche à suivre. Cet article suppose également que vous avez suivi toutes les étapes précédentes, et que vous avez fait fonctionner votre premier run WRF comme cela a été décrit dans l'article 5 - Notre premier run WRF.


Les fichiers GRIB et NetCDF


Nous avons jusqu'ici utilisé des fichiers GRIB provenant du modèle GFS, et après bien du labeur, produit des fichiers de sortie au format NetCDF contenant le résultat de notre simulation météo. Nous ne nous sommes pas trop préoccupés de leur contenu : nous allons maintenant nous familiariser un peu avec ces formats. Pour inspecter un fichier, NCL fournit la commande ncl_filedump qui s'utilise de la façon suivante :

$ ncl_filedump <nom_fichier.grib2>
$ ncl_filedump <nom_fichier.nc>

Essayons par exemple sur l'un de notre fichier GFS. Pour plus de confort, la sortie étant généralement plus longue que l'écran, mieux vaut piper vers la commande less ou un fichier texte :

$ ncl_filedump gfs.t00z.pgrbf00.grib2 | less

Voici un petit extrait de ce que vous obtiendrez :

 Copyright (C) 1995-2013 - All Rights Reserved
 University Corporation for Atmospheric Research
 NCAR Command Language Version 6.1.2
 The use of this software is governed by a License Agreement.
 See http://www.ncl.ucar.edu/ for more details.

Variable: f
Type: file
filename:       gfs.t00z.pgrbf00
path:   gfs.t00z.pgrbf00.grib2
   file global attributes:
   dimensions:
      lat_0 = 181
      lon_0 = 360
      lv_ISBL0 = 26
      lv_AMSL1 = 3
      lv_DBLL2 = 4
      lv_PVL3 = 2
      lv_ISBL4 = 21
      lv_SIGL5 = 4
      lv_ISBL6 = 2
      lv_ISBL7 = 6
   variables:
      float TMP_P0_L1_GLL0 ( lat_0, lon_0 )
         center :       US National Weather Service - NCEP (WMC)
         production_status :    Operational products
         long_name :    Temperature
         units :        K
         _FillValue :   1e+20
         grid_type :    Latitude/longitude
         parameter_discipline_and_category :    Meteorological products, Temperature
         parameter_template_discipline_category_number :        ( 0, 0, 0, 0 )
         level_type :   Ground or water surface
         level :         0
         forecast_time :        0
         forecast_time_units :  hours
         initial_time : 08/15/2007 (00:00)


On retrouve principalement deux sections :
  • la liste des dimensions utilisées pour la grille horizontale (lat_0 et lon_0) et les niveaux verticaux (les variables lv_*).
  • la liste des variables, leur type, les dimensions sur lesquelles elles sont décrites et leur description (unité, etc...) . Des variables du même nom que les dimensions existent permettant de connaître les coordonnées et les niveaux où sont définis chaque point de grille. 
Dans notre exemple, la variable s'appelle TMP_P0_L1_GLL0. Le nommage est un peu barbare, mais c'est parce que le format GRIB ne fournit pas de nom unique pour une variable donnée, NCL applique alors un algorithme pour les nommer de manière unique en fonction de leur type et de leurs dimensions. Voyons plus en détail ce que nous apprend la description :
  • float : notre variable est un nombre de type flottant 32 bits.
  • (lat_0, lon_0) : nous indique que cette variable est à deux dimensions qui sont la latitude et la longitude, ce qui donne une taille de grille de 181x360 selon la définition de lat_0 et lon_0.
  • center et production_status : Différentes informations sur le centre de données qui a créé ce fichier, on sait que c'est un fichier de production opérationnelle du NCEP.
  • Le nom complet de la variable : il s'agit d'une température
  • Units : cette température est en Kelvins.
  • _FillValue : la valeur de remplissage par défaut des données
  • level_type : il s'agit ici d'une variable 2D représentant la surface du sol ou de la mer. On en déduit qu'il s'agit de la température de surface. On peut autrement trouver des variables définies par une hauteur en mètres au-dessus du niveau de la mer, ou encore une altitude isobare.
  • level : altitude du niveau, dans le cas d'une variable 2D uniquement. Par exemple cette information vaudrait 2 pour une température à 2m.
  • forecast_time et forecast_time_units : donne l'heure depuis l'analyse, et dans quelle unité cette heure est exprimée. Ici, il s'agit de 0 puisque c'est l'analyse.
  • initial_time : nous indique la date et l'heure de l'analyse et donc du début de la simulation.
Je vous invite à vous familiariser avec la sortie de cette commande et d'apprendre à repérer les variables que vous avez l'habitude d'analyser (température à 2m, vent à 10m, champ d'humidité 3D etc...). Faites de même avec le fichier wrfout que vous avez obtenu à l'issue de votre simulation. Attention vous devrez le renommer et lui rajouter une extension ".nc" pour que la commande le reconnaisse bien comme un fichier NetCDF. Vous verrez alors que les fichiers NetCDF sont bien plus sympa pour le nommage des variables, par exemple la température à 2 mètres sera décrite comme suit :

      float T2 ( Time, south_north, west_east )
         FieldType :    104
         MemoryOrder :  XY
         description :  TEMP at 2 M
         units :        K
         stagger :     
         coordinates :  XLONG XLAT


Avouez que c'est quand même plus facile de s'y retrouver ! Les informations affichées sont un peu différentes, mais le principe reste le même aussi je ne rentrerais pas dans les détails.

Un premier script


Nous allons maintenant voir quelques bases du langage, et comme nous en sommes à nous familiariser avec les variables nous verrons comment en afficher le contenu. Cela vous sera utile pour débugger, mais aussi pour connaître les hauteurs des niveaux verticaux de votre fichier.  Sans plus attendre, voici le script.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"

begin
    ; Chargement du fichier de sortie WRF
    a = addfile("wrfout.nc","r")  ; Open a file

    ; Récuperation de la liste des heures contenue dans le fichier
    times = wrf_user_getvar(a,"times",-1)

    ; Affichage du contenu
    print(times)
end


Recopiez ce script dans un fichier nommé wrf_premier_script.ncl et tapez ensuite la commande suivante pour l'exécuter (renommez votre fichier de données WRF en wrfout.nc au préalable) :

$ ncl wrf_premier_script.ncl

La sortie de cette commande vous donne la liste des dates contenues dans ce fichier, qui correspondent à notre simulation :

 Copyright (C) 1995-2013 - All Rights Reserved
 University Corporation for Atmospheric Research
 NCAR Command Language Version 6.1.2
 The use of this software is governed by a License Agreement.
 See http://www.ncl.ucar.edu/ for more details.


Variable: times
Type: string
Total Size: 200 bytes
            25 values
Number of Dimensions: 1
Dimensions and sizes:    [25]
Coordinates:
Number Of Attributes: 2
  _FillValue :    missing
  description :    times in file
(0)    2013-07-28_00:00:00
(1)    2013-07-28_01:00:00
(2)    2013-07-28_02:00:00
(3)    2013-07-28_03:00:00
(4)    2013-07-28_04:00:00
(5)    2013-07-28_05:00:00
(6)    2013-07-28_06:00:00
(7)    2013-07-28_07:00:00
(8)    2013-07-28_08:00:00
(9)    2013-07-28_09:00:00
(10)    2013-07-28_10:00:00
(11)    2013-07-28_11:00:00
(12)    2013-07-28_12:00:00
(13)    2013-07-28_13:00:00
(14)    2013-07-28_14:00:00
(15)    2013-07-28_15:00:00
(16)    2013-07-28_16:00:00
(17)    2013-07-28_17:00:00
(18)    2013-07-28_18:00:00
(19)    2013-07-28_19:00:00
(20)    2013-07-28_20:00:00
(21)    2013-07-28_21:00:00
(22)    2013-07-28_22:00:00
(23)    2013-07-28_23:00:00
(24)    2013-07-29_00:00:00


Voyons maintenant les explications de rigueur :
  • un certain nombre de lignes "load" chargent des librairies contenant les définitions de toutes les fonctions que nous utiliserons pour nos cartes et graphiques. Cette liste est à peu près exhaustive et vous donnera accès à tout le panel nécessaire pour toute utilisation.
  • le corps du script est délimité par begin et end.
  • une ligne = une instruction, si toutefois votre instruction est trop longue vous pouvez la répartir sur plusieurs lignes en mettant un antislash en fin de ligne.
  • les commentaires sont marqués par le point virgule
  • les variables sont déclarées au moment de leur affectation
  • la fonction addfile ouvre un fichier, ici en lecture
  • la fonction wrf_user_getvar est une fonction utilitaire permettant d'obtenir le contenu d'une variable du fichier pour une échéance donnée ou pour toutes les échéances (-1). Elle permet également de renvoyer des variables diagnostic qui ne font pas partie du fichier et sont recalculées (theta, etc...). Cette fonction fait partie de toute une bibliothèque permettant de simplifier l'utilisation des données issues de WRF ARW (et uniquement ce modèle - nous le verrons quand nous aborderons NMM) en gérant par exemple de manière transparente le type de grille un peu spécifique de WRF. 
  • La fonction print permet d'afficher le contenu d'une variable, ou un message. Ici, il s'agit de la liste des heures de prévision, lue précédemment dans notre fichiers.

Notre première carte


Mais enchainons sans plus attendre sur un script au résultat plus visuel ! Nous allons créer un script qui affiche les isobares de pression au niveau de la mer. Enregistrez le code source ci-dessous dans un fichier NCL et exécutez-le.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"

begin
    ; Chargement du fichier de sortie WRF
    a = addfile("wrfout.nc","r")
    times = wrf_user_getvar(a,"times",-1)

    ; Ouverture du workspace graphique
    type = "x11"
    wks = gsn_open_wks(type,"pression")

    ; Ressources communes aux différentes couches
    res = True
    res@MainTitle = "METEO BLOIS"

    ; Parcours des previsions
    ntimes = dimsizes(times)
    do it = 0,ntimes-1,1
   
        res@TimeLabel = times(it)

        ; Tracé de la pression de surface   
        press = wrf_user_getvar(a,"slp", it)

        contour = wrf_contour(a,wks,press,res)

        ; Dessin de la carte en overlay
        pltres = True
 

        ; On précise la couleur des traits de la carte pour plus de lisibilité
        mpres = True
        mpres@mpGeophysicalLineColor      = "Black"
        mpres@mpGridLineColor             = "Black"
        mpres@mpLimbLineColor             = "Black"
        mpres@mpNationalLineColor         = "Black"
        mpres@mpPerimLineColor            = "Black"

        
        ; On ajoute notre plot de pression sur la carte
        plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
  
    end do
end


Le programme NCL ouvre alors une fenêtre qui affiche la carte comme illustré ci-dessous. Cliquez dans la fenêtre pour dérouler les échéances.

Contours de pression au niveau de la mer
Contour de pression au niveau de la mer

La fonction gsn_open_wks permet d'ouvrir un workspace graphique, c'est-à-dire l'objet dans lequel on va dessiner. Cet objet peut être de plusieurs types, essentiellement X11 (affichage écran) et PS (postscript, pour conversion en image). On lui donne un nom de fichier qui servira pour les sorties postscript.

Nous définissons ensuite le titre de notre carte. Cela se fait grâce à une variable NCL que nous décorons d'un certain nombre de propriétés appelées ressources. Les ressources sont accessibles avec la syntaxe variable@ressource. Les différentes fonctions graphiques regorgent de ressources permettant de paramétrer l'aspect visuel souhaité : nous n'en verrons dans cet article qu'un très petit nombre. Dans cet exemple, nous laissons le plus de paramètres d'origine que possible pour simplifier, seules les couleurs de la carte ont été forcées pour plus de lisibilité.

Le parcours de notre fichier se fait par la boucle "do ... end do", équivalent de la boucle "for" de bien des langages. On retrouve les trois parties classiques, à savoir l'initialisation de la variable de boucle, la condition de bouclage et l'incrément de la variable de boucle. La fonction dimsizes donne la taille de notre tableau nous permettant de borner notre boucle. Une seule valeur est ici renvoyée puisque notre tableau des heures n'a qu'une dimension.

On retrouve la désormais connue fonction wrf_get_var pour obtenir la variable pression au niveau de la mer à l'échéance souhaitée. Le tracé des isobares se fait avec la fonction wrf_contour, qui prend en paramètre notre objet fichier, l'objet workspace, la variable à tracer et les ressources de paramétrage visuel.

Enfin, nous utilisons la fonction wrf_map_overlays pour tracer la carte avec notre contour en superposition. Celle-ci prend un tableau d'objets graphiques à superposer : ceci est marqué par la syntaxe "(/contour/)". Il est ainsi possible de définir autant de couches que nous le souhaitons. Cette fonction finalise également le tracé du graphique, en lançant le dessin effectif de nos objets sur l'écran (ou le fichier) et passe à la frame suivante. Lors d'un affichage écran, NCL attend que l'utilisateur clique pour passer à la frame suivante.



Ajout du géopotentiel 500 hPa et améliorations


Par défaut, NCL choisit les contours de manière à obtenir un espacement homogène. Dans notre exemple, cela résulte en un espacement de 2 hPa entre les isobares, ce qui est assez bien mais pas forcément ce que l'on souhaite : nous allons les espacer de 1 hPa. Pour cela, nous utilisons les ressources suivantes :

optpres@cnLevelSelectionMode = "ManualLevels"
optpres@ContourParameters    =  1


La ressource cnLevelSelectionMode spécifie que nous souhaitons définir manuellement les niveaux de notre tracé de contour, et la ressource ContourParameters est un raccourci proposé par la fonction wrf_contour permettant de définir l'espacement entre les isolignes : ici il s'agit d'un espacement de 1 unité.

Afin d'agrémenter encore un peu notre carte, nous allons rajouter une analyse d'altitude en superposant le géopotentiel 500 hPa à notre analyse de surface. Cette donnée n'est pas disponible directement dans notre fichier, nous devrons donc l'interpoler à partir des données présentes : la pression sur les niveaux du modèle et l'altitude (z) des niveaux :

z = wrf_user_getvar(a,"z", it)
p = wrf_user_getvar(a, "pressure", it)


Pour interpoler, nous utilisons la fonction wrf_user_intrp3d qui réalise une coupe des données d'une grille selon un plan horizontal ou vertical donné.

g500 = wrf_user_intrp3d(z, p, "h", 500, 0., False)

Voici la signification des paramètres :
function wrf_user_intrp3d (
 var3d      : numeric,  
 vert       : numeric,  
 plot_type  : string,   
 loc    [*] : numeric,  ; up to four values
 angle      : numeric,  
 res    [1] : logical   
)
  • var3d = z : la variable 3D que l'on veut interpoler. On veut en effet obtenir la hauteur à pression constante.
  • vert = p : la coordonnée verticale à utiliser. Cela peut être une hauteur ou une pression, dans notre cas on travaille en coordonnée pression.
  • plot_type = "h" : indique le type de coupe que l'on veut obtenir, ici il s'agit d'une coupe horizontale.
  • loc = 500 : paramètre de la coupe, une seule valeur est nécessaire ici pour notre coupe horizontale qui correspond à nos 500 hPa.
  • angle = 0 : indique l'angle de la coupe, valable seulement pour les coupes verticales. On laisse ici 0.
  • res = False : sert à indiquer qu'une coupe s'effectue entre un point A et un point B, ce qui n'est pas le cas ici.
On rajoute ensuite une description de la variable ainsi créée pour la légende :

g500@description = "Geopotentiel 500 hPa"

Le tracé se réalise toujours avec la fonction wrf_contour, mais cette fois-ci nous paramétrons les contours avec un remplissage de couleur. L'espacement est également précisé ainsi que les valeurs mini et maxi de tracé, ce qui permet de "figer" les couleurs par rapport à la palette et de ne pas avoir une échelle qui varie à chaque carte en fonction des valeurs fournies par le modèle.

optg500 = res
optg500@cnFillOn              = True
optg500@gsnSpreadColors     = True
optg500@cnLevelSelectionMode = "ManualLevels"
optg500@cnMinLevelValF     = 5000.
optg500@cnMaxLevelValF     = 6040.
optg500@cnLevelSpacingF     = 40.
optg500@ContourParameters    = 40

geop = wrf_contour(a, wks, g500, optg500)


Il ne reste plus qu'à effectuer l'overlay des deux plots :

plot = wrf_map_overlays(a,wks,(/geop, contour/),pltres,mpres)

Voici le code source complet et le résultat attendu.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"

begin
    ; Chargement du fichier de sortie WRF
    a = addfile("wrfout.nc","r")
    times = wrf_user_getvar(a,"times",-1)

    ; Ouverture du workspace graphique
    type = "x11"
    wks = gsn_open_wks(type,"plt_press_geop500")

    ; Ressources communes aux différentes couches
    res = True
    res@MainTitle = "METEO BLOIS"
    res@Footer = False

    ; Parcours des previsions
    ntimes = dimsizes(times)
    do it = 0,ntimes-1,1
   
        res@TimeLabel = times(it)

        ; Tracé de la pression de surface   
        press = wrf_user_getvar(a,"slp", it)
        press@description     = "Pression au niveau de la mer"

        optpres = res
        optpres@cnFillOn              = False
        optpres@cnLevelSelectionMode = "ManualLevels"
        optpres@ContourParameters    =  1

        contour = wrf_contour(a,wks,press,optpres)

        ; Tracé du géopotenciel 500 hpa
        z = wrf_user_getvar(a,"z", it)
        p = wrf_user_getvar(a, "pressure", it)
        g500 = wrf_user_intrp3d(z, p, "h", 500, 0., False)
        g500@description = "Geopotentiel 500 hPa"

        optg500 = res
        optg500@cnFillOn              = True
        optg500@gsnSpreadColors     = True
        optg500@cnLevelSelectionMode = "ManualLevels"
        optg500@cnMinLevelValF     = 5000.
        optg500@cnMaxLevelValF     = 6040.
        optg500@cnLevelSpacingF     = 40.
        optg500@ContourParameters    = 40

        geop = wrf_contour(a, wks, g500, optg500)

        ; Dessin de la carte en overlay
        pltres = True
        pltres@LatLonOverlay = False

        mpres = True
        mpres@mpGeophysicalLineColor      = "Black"
        mpres@mpGridLineColor             = "Black"
        mpres@mpLimbLineColor             = "Black"
        mpres@mpNationalLineColor         = "Black"
        mpres@mpPerimLineColor            = "Black"

        plot = wrf_map_overlays(a,wks,(/geop, contour/),pltres,mpres)      
    end do
end


WRF Pression et géopotentiel à 500 hPa.
Pression et géopotentiel à 500 hPa.



Transformer en image


Bien souvent, on souhaite utiliser nos données pour obtenir des images, pour les afficher sur un site web par exemple. Pour cela, nous allons tout d'abord effectuer le rendu en postscript en modifiant dans notre script le type de workspace de X11 en PS.

; Ouverture du workspace graphique
type = "ps"
wks = gsn_open_wks(type,"plt_press_geop500")


NCL génèrera alors un fichier PostScrip nommé "plt_press_geop500.ps" contenant une page pour chaque heure du fichier WRF. Ce fichier devra être converti en une série d'images PNG à l'aide de l'utilitaire GhostScript :

$ gs -sDEVICE=png16m -dNOPAUSE -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4 -r100 -sOutputFile="plt_press_geop500_%03d.png" -c "<</Orientation 0>> setpagedevice" --f plt_press_geop500.ps quit

La commande est un peu complexe mais en voici les paramètres essentiels :
  • -sDEVICE=png16m : choix du format de sortie, du PNG en 16 millions de couleurs.
  • -sOutputFile="plt_press_geop500_%03d.png" : nommage de notre fichier de sortie. Le "%03d" est un marqueur indiquant à Ghostscript qu'il devra numéroter les fichiers de sortie.
  • -c "<</Orientation 0>> setpagedevice" : pas de rotation de la page. Peut être utile si votre graphique est rendu en format paysage pour le remettre dans le bon sens.
  • -r100 : on veut faire un rendu avec une résolution de 100dpi, ce qui est approximativement voire un peu mieux que la définition d'un écran.

 

Conclusion


Il y aurait encore bien des choses à dire sur NCL tant ce système est riche, nous n'avons fait qu'en survoler les possibilités. Mais ce n'est pas l'objectif premier de cet article, et la place me manquerait pour tout détailler. Pour approfondir le lecteur est invité à se référer à l'abondante documentation du langage NCL, accompagnée de nombreux exemples. Après cette courte pause "ludique", le prochain article de cette série de tutoriaux sera dédié cette fois à la mise en oeuvre de WRF NMM : Partie 7 - Mise en oeuvre de WRF-NMM.