Python wrfout文件提取气象站点数据

一、数据读取

import numpy as np
import datetime
from wrf import getvar, ALL_TIMES, ll_to_xy, interplevel, to_np
import matplotlib.pyplot as plt
import netCDF4 as nc

wrfout_dir = "./data/wrfout_d02_2017-08-22_12"
wrfout     = nc.Dataset(wrfout_dir, mode="r")

1.1 查看和获取变量

获取所有变量名

wrfout.variables.keys()

获取单个数据

wrfout.variables["T2"]

wrfout['T2']

1.2 查看和获取维度

wrfout.dimensions

wrfout.dimensions['south_north']

wrfout.dimensions['south_north'].size

1.3 质量点和跳点mass/stag

wrfout.dimensions['south_north'].size

wrfout.dimensions['south_north_stag'].size

1.4 getvar函数

lat2D      = getvar(wrfout, "lat",     )# units: decimal degrees
lon2D      = getvar(wrfout, "lon",   meta=False  )  # units: decimal degrees
times      = getvar(wrfout, "times", meta=False, timeidx=ALL_TIMES) 

#读取3维u，v， 位势高度和地形高度，特殊属性用getvar来读取
ua     = getvar(wrfout, "ua" , timeidx=ALL_TIMES, meta=False, units="m s-1")
va     = getvar(wrfout, "va" , timeidx=ALL_TIMES, meta=False, units="m s-1")
height = getvar(wrfout, "z"  , timeidx=ALL_TIMES, meta=False, units="m"    )
ter1   = getvar(wrfout, "ter", timeidx=0,         meta=False, units="m"    )
ter2   = getvar(wrfout, "ter", timeidx=0,         meta=True , units="m"    )

# U10 这种本来就有的属性，可以直接按照nc文件的方式来读取
u10 = wrfout.variables['U10'][:]
v10 = wrfout.variables['V10'][:]
wspd_10m = [] 

---

二、站点数据提取

lonSta1 = 114.1  # TC中心站点经纬度
latSta1 = 21.8

lonSta2 = 116.0 # 外围站点经纬度
latSta2 = 22.0

#返回的iSta, jSta为X，Y方向的索引
iSta1, jSta1 = ll_to_xy(wrfout, latSta1, lonSta1, timeidx=0, meta=False)
iSta2, jSta2 = ll_to_xy(wrfout, latSta2, lonSta2, timeidx=0, meta=False)
iSta1, jSta1

#10m风速的提取
wspd_10m_sta1 = [] 
wspd_10m_sta2 = [] 
for it in range(len(times)):
    u10_it = u10[it, jSta1, iSta1]
    v10_it = v10[it, jSta1, iSta1]
    wspd_10m_sta1.append(np.sqrt(u10_it**2 + v10_it**2))
    
    u10_it = u10[it, jSta2, iSta2]
    v10_it = v10[it, jSta2, iSta2]
    wspd_10m_sta2.append(np.sqrt(u10_it**2 + v10_it**2))

---

三、数据可视化

# wrf time转为北京时
def wrftime_to_CST_date(wrftime):
    datestr_wrf  = str(wrftime)[0:13] # 2017-08-22T12, UTC
    datetime_wrf = datetime.datetime.strptime(datestr_wrf, '%Y-%m-%dT%H')
    datetime_CST = datetime_wrf + datetime.timedelta(hours=8)
    datestr_CST  = datetime_CST.strftime('%Y-%m-%d_%H')
    return datestr_CST

x = range(len(times))
plt.figure(figsize=(8,6),dpi=100)
plt.plot(x, wspd_10m_sta1,    linewidth=2.0, color='red',   label='sta1'   )
plt.plot(x, wspd_10m_sta2,    linewidth=2.0, color='blue',  label='sta2'   )

plt.legend(fontsize=20)

init_date_cst = wrftime_to_CST_date(times[0])
title_name = ' init_time: '+ init_date_cst 
plt.title(title_name, fontsize=20)
plt.xlabel('time(hour)', fontsize=20)
plt.ylabel('wind speed(m/s)', fontsize=20)
fig_name = 'WRF_station_wind_'+init_date_cst+'.png'
plt.savefig(fig_name)