用Pygrib读取grb2文件，读到的值跟wgrib2的不同

不远行

处理的是中尺度数值天气预报的的grb2文件，是由GFS预测资料降尺度的来的。现在想从中提取出想要的风速、2m温度等参量制作天气预报，

开始用pygrib定义对象，提取grb2文件中特定坐标处的各个参量，后来尝试用wgrib2工具读取相同的文件，发现读到的数值不一样。下面是python脚本：

#!/usr/bin/env python
#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.use('AGG')
from pylab import *
import numpy as np
#import matplotlib
import Image
import pylab as m
import pygrib
import commands
#import pygrib
#rc('font',**{'family':'serif','serif':['simhei']})


# COMPROBAR LOS ARGUMENTOS DE ENTRADA
if int(size(sys.argv)) < 4:
        print "Faltan argumentos:"
        print "./GRIBtoAscii.py y x File"
        exit()


j = int(sys.argv[1])
i = int(sys.argv[2])
Fichero=sys.argv[3]


# open a GRIB2 file, create a Grib2 class instance.
grbs = pygrib.open(str(Fichero))  
grbs.seek(0)


#LON = grbs[92].values[25,:]-360
#LAT = grbs[91].values[:,0]


##FUNCION ENCONTRAR LAT/LON + CERCANA
#def find_nearest(array,value):
    #idx=(np.abs(array-value)).argmin()
    #return array[idx], idx


##ENCONTRAR I,J DE LA LAT/LON
#NearLon=find_nearest(LON,longitude)
#NearLat=find_nearest(LAT,latitude)


#SACAR VARIABLES DEL FICHERO GRIB Y CALCULAR FECHA DE LA PREDICCION
StartDate = grbs[1].dataDate
Alc = grbs[1].endStep
ForecastDate = commands.getoutput('date -u -d "'+str(Alc)+' minutes "'+str(StartDate)+'" 20 UTC" "+%Y%m%d%H"')

P = grbs[6].values[j,i]
T2M = grbs[3].values[j,i]


U10M = grbs[1].values[j,i]
V10M = grbs[2].values[j,i]


RH2M = grbs[5].values[j,i]
RAIN = grbs[7].values[j,i]


CLOUD = grbs[19].values[j,i]


#MODULOS DE LA VELOCIDAD
Vel10M = sqrt(U10M*U10M+V10M*V10M)       


#DIRECCION DEL VIENTO PREDICHO
def CalcDirection(U,V):


        if U==0. and V==0.:
                Dir=0.
        else:
                alfa=math.atan2(float(abs(U)),float(abs(V)))*(180/pi)


        if U>0. and V>0.:
                Dir=180+alfa
               
        elif U>0. and V<0.:
                Dir=360-alfa
               
        elif U<0. and V>0.:
                Dir=180.-alfa
               
        else:
                Dir=alfa
        return Dir


Dir10M = CalcDirection(U10M,V10M)




print "%s\t%i\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f" \
% (str(ForecastDate), Alc, T2M-273.15, U10M, V10M, Vel10M, Dir10M, P/100, RH2M, RAIN, CLOUD)


grbs.close()


wgrib2这个工具目前不太会用，用的是简单指令： wgrib2 xxx.grb2 -s -ijlat 225 435


不知道为什么两种办法读出来的值不一样。

阿阿飞飞

pygrib使用的是ECMWF的grib API接口，不知道是不是跟GFS的不太一致。
一般来说wgrib给出的结果都是正确的。你在linux系统下，可以使用很多其它函数，如pynio、iris等

delete99999

我还不会这个，下载下来研究一下

ntwd87

楼主用的cygwin？

dolphin0520

谢谢分享！

余师兄

subtropical

楼主，请问您知道怎么用pygrib把模式资料grb转为文本文件吗