Discussion :

[flyingfr]

Bonjour,

Je dispose d'un fichier grib avec des valeurs de vent fonction de la latitude, de la longitude et du temps (0,+3,......).

J'ouvre sans probleme ce fichier avec

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
Fichier_grib='Gribs/test_simple.grb'
grbs= pygrib.open(Fichier_grib)

je peux imprimer la liste des messages avec

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
grbs.seek(0)
print (' resultat de grbs.seek(0)')
print()
for grb in grbs:
    print(grb)

et j'ai bien dans les réponses par exemple :

.....
18:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10:fcst time 0 hrs:from 201702180000
19:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10:fcst time 3 hrs:from 201702180000
.....
36:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10:fcst time 3 hrs:from 201702180000
....


je voudrais récupérer les valeurs '10 metre U wind component' et '10 metre V wind component' par exemple pour le temps +3h la latitude -32 et la longitude 55.

Comment puis-je faire?

Merci par avance

[marco056]

Pour t'aider, il faudrait peut-être le fichier. Tu télécharges où ?

[flyingfr]

Bonjour,

A force de tâtonnements, voilà la solution pour le prochain qui aurait ce genre de question !
On récupère un tableau indexé à trois dimensions sur le paramètre cherché à partir du fichier par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
fileidx=pygrib.index(Fichier_grib,'name')
u= fileidx.select(name="10 metre U wind component")

Il ne reste plus qu'à extraire les valeurs avec les indices entiers commençant à 0 0 pour le coin en haut à gauche avec un balayage horizontal suivant le pas du grib ( ici pas de 0.25° )
l'indice de temps correspond à l'indice de l'heure ( grib toutes les 3 h pour moi )

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

 vu0=u[itemps0]['values'][ilat][jlong]

Ce qui donne en global avec une interpolation entre 2 indices de temps

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Feb 19 09:00:37 2017
 
@author: jph
"""
 
import pygrib
import numpy as np
import datetime
import time
import math
 
 
Fichier_grib='Gribs/position_19_02.grb'
grbs= pygrib.open(Fichier_grib)
fileidx=pygrib.index(Fichier_grib,'name')
print()
print (' Nom du Fichier Grib utilisé  : ', grbs.name)
 
 
g=grbs[1]
 
day     =g['day']
month   =g['month']
year    =2000+g['yearOfCentury']
hour    =g['hour']
minute  =g['minute']
second =g['second']
 
latfirst=g['latitudeOfFirstGridPointInDegrees']
lonfirst=g['longitudeOfFirstGridPointInDegrees']
latlast =g['latitudeOfLastGridPointInDegrees']
lonlast =g['longitudeOfLastGridPointInDegrees']
 
paslat=g['iDirectionIncrementInDegrees']
paslong=g['jDirectionIncrementInDegrees']
 
print(" Heure et date du grib   {}/{}/{} -  {}h {}mn {}s ".format(day,month,year,hour,minute,second ))
print (" Grib entre longitudes de {} à {} et latitudes de {} à {}  ".format(lonfirst,lonlast,latfirst,latlast))
print (" Pas du Grib en latitude : {}° longitude : {}° ".format(paslat,paslong))
print (' Nombre de messages  : ', grbs.messages)
print ()
 
 
 
# Date et lieu choisis pour la prévision
lon_p=28.5
lat_p=-57.75
year_p=2017
month_p=2
day_p=20
hour_p=21
mn_p=0
 
date_p=time.mktime((year_p,month_p,day_p,hour_p,mn_p,0,0,0,0))  # Instant de la prevision en s 
date_i=time.mktime((year,month,day,hour,minute,0,0,0,0))        # Instant 0 du Grib
 
 
 
def fvent(longitude,latitude,date_p):
 
 
    ecart=date_p-date_i   # ecart avec l'instant 0 du grib 
 
    itemps0=int(ecart/ 3600//3) #indice du temps arrondi inferieur
    itemps1=(ecart/ 3600/3)    #indice du temps non arrondi
 
    jlong=int(abs(longitude-lonfirst)/paslat) #indice de la longitude
    ilat=int(abs( latitude-latfirst)/paslat) #indice de la latitude
 
    u= fileidx.select(name="10 metre U wind component")
    v=fileidx.select(name="10 metre V wind component")
 
    vu0=u[itemps0]['values'][ilat][jlong]
    vu1=u[itemps0+1]['values'][ilat][jlong]   
    vv0=v[itemps0]['values'][ilat][jlong]
    vv1=v[itemps0+1]['values'][ilat][jlong]
 
    # interpolation
    vu=vu0 +(vu1-vu0)*(itemps1-itemps0)
    vv=vv0 +(vv1-vv0)*(itemps1-itemps0)
 
    vn=math.sqrt(vu**2+vv**2)*1.94384# Vitesse du vent en noeuds
    av=270-math.atan(vv/vu)*180/3.1416# Angle du vent en degrés
    date=time.strftime('%Hh%M le %d/%m/%Y',time.localtime(date_p)) # Temps formaté 
    return (longitude,latitude,date,vn,av)
 
 
 
 
vent =fvent(lon_p,lat_p,date_p)
print()
 
 
print (" En {}° de longitude et {}° de latitude à {} Vent de :{:6.2f} Noeuds  Direction {:6.2f}° ".format(vent[0],vent[1],vent[2],vent[3],vent[4]))

Maintenant je vais m'attaquer à la représentation sous matplotlib.

A+

[flyingfr]

Bonjour,

Pour la question de marco56 sur comment je récupère mes fichiers Gribs :

J''utilise Zygrib qui récupère pour moi le fichier auprès de la NOAA correspondant directement à la zone et aux paramètres souhaités.
Je pourrais aussi le récupérer directement par mail, mais il faut définir les paramètres, c'est plus facile avec Zygrib.

Cordialement.
A+

[marco056]

Ok