python3读取fnl 1*1资料双线性水平插值到格点，垂直对数插值到43层

Lucian

python3在ubuntu中读取fnl 1*1资料中的温度和相对湿度，把相对湿度转化成ppmv。双线性水平插值到指定经纬度，垂直对数插值到rttov9.3计算radiance的43层

1. import numpy as np
   
2. from scipy import interpolate
   
3. import pygrib
   
4. 
   
5. def rh2q_ppmv(RH,P,T):
   
6.     '''Conversion reletive humidity to lnq in ppmv
   
7.        RH-%, P-hPa, T-K'''
   
8.     t = T-273.15
   
9.     if T <=233.15:
   
10.         e_s=6.1078*np.exp(21.8746*(t)/(T-7.66))
    
11.     elif T>=258.15:
    
12.         e_s=6.1078*np.exp(17.269*(t)/(T-35.86))
    
13.     else:
    
14.         e_s_water=e_s=6.1078*np.exp(17.269*(t)/(T-35.86))
    
15.         e_s_ice=6.1078*np.exp(21.8746*(t)/(T-7.66))
    
16.         e_s=0.002*((80+2*t)*e_s_water-(30+2*t)*e_s_ice)
    
17.    
    
18.     e = RH/100*e_s # hPa   
    
19.     q_ppmv = 1000000*e/P
    
20.     return q_ppmv
    
21. 
    
22. # 指定fnl资料路径
    
23. grbs = pygrib.open('/home/xxxxx/Documents/fnl/fnl_20180110_00_00.grib2')
    
24. grbs.seek(0)
    
25. 
    
26. #for key in grbs:
    
27.     #print(key)
    
28. 
    
29. Temperatures = grbs.select(name='Temperature')[:] # 读取31层温度
    
30. RHs = grbs.select(name='Relative humidity')[:] # 读取31层相对湿度
    
31. 
    
32. lats, lons = Temperatures[0].latlons() # 读取经纬度
    
33. lats, lons = lats[:,0], lons[0,:] # 压缩经纬度到一维
    
34. 
    
35. # fnl资料里31层气压值
    
36. fnl_P_levels_hpa = [1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 20.0, 30.0, 50.0, 70.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 925.0, 950.0, 975.0, 1000.0]
    
37. fnl_lnP_levels_hpa = np.log(fnl_P_levels_hpa) # 取对数
    
38. 
    
39. fnl_Vertical_T = []
    
40. fnl_Vertical_RH = []
    
41. fnl_Vertical_q_ppmv = []
    
42. for key in range(31):
    
43.     fx = interpolate.interp2d(lons, lats, Temperatures[key].values, kind='linear')
    
44.     T = fx(120.5,25.5)# 将fnl31层温度插值到120.5经度，北纬25.5度   
    
45.     fx = interpolate.interp2d(lons, lats, RHs[key].values, kind='linear')
    
46.     RH = fx(120.5,25.5)# 将fnl31层相对湿度插值到120.5经度，北纬25.5度
    
47. 
    
48.     q_ppmv_cal = rh2q_ppmv(RH,fnl_P_levels_hpa[key],T)
    
49.     fnl_Vertical_T.append(float(T))
    
50.     fnl_Vertical_RH.append(float(RH))
    
51.     fnl_Vertical_q_ppmv.append(float(q_ppmv_cal))
    
52.    
    
53. #print(fnl_Vertical_q_ppmv)
    
54. 
    
55. rttov93_P_levels = [0.1, 0.29, 0.69, 1.42, 2.611, 4.407, 6.95, 10.37, 14.81, 20.4, 27.26, 35.51, 45.29, 56.73, 69.97, 85.18, 102.05, 122.04, 143.84, 167.95, 194.36, 222.94, 253.71, 286.6, 321.5, 358.28, 396.81, 436.95, 478.54, 521.46, 565.54, 610.6, 656.43, 702.73, 749.12, 795.09, 839.95, 882.8, 922.46, 957.44, 985.88, 1005.43, 1013.25]
    
56. rttov93_lnP_levels = np.log(rttov93_P_levels) # 取对数
    
57. 
    
58. fx = interpolate.interp1d(fnl_lnP_levels_hpa, fnl_Vertical_T, kind='linear', fill_value='extrapolate' )
    
59. rttov93_Vertical_T = fx(rttov93_lnP_levels) # 温度插值
    
60. 
    
61. fx = interpolate.interp1d(fnl_lnP_levels_hpa, fnl_Vertical_q_ppmv, kind='linear', fill_value='extrapolate' )
    
62. rttov93_Vertical_q_ppmv = fx(rttov93_lnP_levels) # q_ppmv插值
    
63. rttov93_Vertical_lnq_output = np.log(rttov93_Vertical_q_ppmv/1000) # 转化成rttov93的存储格式
    
64. # 插值结果
    
65. print(rttov93_Vertical_T)
    
66. print(rttov93_Vertical_lnq_output)

复制代码

Lucian

二楼占楼自用

yuchen

学习了，感谢楼主！

kust

{:5_275:}zan

独孤酒见

给楼主点赞，分享要大赞！

x6xin

{:5_231:}{:5_231:}{:5_231:}{:5_231:}{:5_231:}

syl136

python可以读grib数据啊？？！！！学习学习！！一直以为不能读grib

rosecjm

感谢分享！收藏，学习~

python游

pygrib 在windows下安装不上，请各位指教

Lucian

python游 发表于 2018-8-18 21:54
pygrib 在windows下安装不上，请各位指教

老哥我的贴子是ubuntu，linux的一种，windows安装pygrib百度有教程的