Python气象数据处理与绘图：WRF模式模拟数据后处理之计算篇

1、导入模块

1import numpy
2import cartopy
3import matplotlib
4from netCDF4 import Dataset
5from xarray import DataArray
6from wrf import getvar, interplevel, vertcross,vinterp, ALL_TIMES, CoordPair, xy_to_ll, ll_to_xy 
7import os
8import warnings
9warnings.filterwarnings('ignore')

2、读取数据（单个和多个）

 1WRF_DIRECTORY = "../input/wrf3880"
 2WRF_FILES = ["wrfout_d01_2005-08-28_00_00_00",
 3             "wrfout_d01_2005-08-28_12_00_00",
 4             "wrfout_d01_2005-08-29_00_00_00"]
 5
 6_WRF_FILES = [os.path.abspath(os.path.join(WRF_DIRECTORY, f)) for f in WRF_FILES]
 7
 8for f in _WRF_FILES:
 9    if not os.path.exists(f):
10        raise ValueError("{} does not exist. "
11            "Check for typos or incorrect directory.".format(f))
12
13def single_wrf_file():
14    global _WRF_FILES
15    return _WRF_FILES[0]
16
17def multiple_wrf_files():
18    global _WRF_FILES
19    return _WRF_FILES

3、读取位势高度（HGT）

1file_path = single_wrf_file()
2wrf_file = Dataset(file_path)
3hgt = getvar(wrf_file, "HGT", timeidx=0)
4hgt

4、读取海平面气压（SLP）

1file_path = single_wrf_file()
2wrf_file = Dataset(file_path)
3slp = getvar(wrf_file, "slp", timeidx=0, units="hPa")
4slp

5、合并所有时刻的海平面气压（cat方法）

1file_paths = multiple_wrf_files()
2wrf_files = [Dataset(f) for f in file_paths]
3slp = getvar(wrf_files, "slp", timeidx=ALL_TIMES, method="cat")
4slp

6、合并所有时刻的海平面气压（join方法）

1file_paths = multiple_wrf_files()
2wrf_files = [Dataset(f) for f in file_paths]
3slp = getvar(wrf_files, "slp", timeidx=ALL_TIMES, method="join")
4slp

7、气压插值到500hPa

1file_path = single_wrf_file()
2wrf_file = Dataset(file_path)
3pres = getvar(wrf_file, "pressure", timeidx=0)
4ht = getvar(wrf_file, "z", timeidx=0, units="dm")
5ht

1ht_500 = interplevel(ht, pres, 500.0)
2ht_500

8、插值到垂直截面（vertcross）

1file_path = single_wrf_file()
2wrf_file = Dataset(file_path)
3bottom_left = CoordPair(x=0, y=0)
4top_right = CoordPair(x=-1, y=-1)
5wspd_wdir = getvar(wrf_file, "wspd_wdir", timeidx=0, units="kt")            
6wspd = wspd_wdir[0,:]
7ht = getvar(wrf_file, "z", timeidx=0)
8wspd_cross = vertcross(wspd, ht, start_point=bottom_left, end_point=top_right, latlon=True)
9wspd_cross

9、插值到Theta-e坐标系（vinterp）

 1file_path = single_wrf_file()
 2wrf_file = Dataset(file_path) 
 3pres = getvar(wrf_file, "pressure", timeidx=0)
 4interp_levels = [280., 285., 290., 292., 294., 296., 298., 300., 305., 310.]
 5pres_eth = vinterp(wrf_file, 
 6                   field=pres, 
 7                   vert_coord="theta-e", 
 8                   interp_levels=interp_levels, 
 9                   extrapolate=True, 
10                   field_type="pressure", 
11                   log_p=False,
12                   timeidx=0)
13pres_eth

10、经纬度和索引的转换（ xy_to_ll, ll_to_xy ）

1file_path = single_wrf_file()
2wrf_file = Dataset(file_path)
3lat_lon = xy_to_ll(wrf_file, [20, 30], [50,75])
4lat_lon

1x_y = ll_to_xy(wrf_file, lat_lon[0,:], lat_lon[1,:])
2print(x_y)