import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import cartopy.io.shapereader as shpreader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from netCDF4 import Dataset
from wrf import to_np, getvar, interplevel, ALL_TIMES, latlon_coords, get_cartopy
import warnings
warnings.filterwarnings('ignore')
3. 准备工作: 兰伯特刻度线补丁子程序
from copy import copy
import numpy as np
import shapely.geometry as sgeom
import cartopy.crs as ccrs
deffind_side(ls, side):"""
Given a shapely LineString which is assumed to be rectangular, return the
line corresponding to a given side of the rectangle.
"""
minx, miny, maxx, maxy = ls.bounds
points ={'left':[(minx, miny),(minx, maxy)],'right':[(maxx, miny),(maxx, maxy)],'bottom':[(minx, miny),(maxx, miny)],'top':[(minx, maxy),(maxx, maxy)],}return sgeom.LineString(points[side])deflambert_xticks(ax, ticks):"""Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
te =lambda xy: xy[0]
lc =lambda t, n, b: np.vstack((np.zeros(n)+ t, np.linspace(b[2], b[3], n))).T
xticks, xticklabels = _lambert_ticks(ax, ticks,'bottom', lc, te)
ax.xaxis.tick_bottom()
ax.set_xticks(xticks)
ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick)for xtick in xticklabels])deflambert_yticks(ax, ticks):"""Draw ticks on the left y-axis of a Lamber Conformal projection."""
te =lambda xy: xy[1]
lc =lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n)+ t)).T
yticks, yticklabels = _lambert_ticks(ax, ticks,'left', lc, te)
ax.yaxis.tick_left()
ax.set_yticks(yticks)
ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick)for ytick in yticklabels])def_lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):"""Get the tick locations and labels for an axis of a Lambert Conformal projection."""
outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
axis = find_side(outline_patch, tick_location)
n_steps =30
extent = ax.get_extent(ccrs.PlateCarree())
_ticks =[]for t in ticks:
xy = line_constructor(t, n_steps, extent)
proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:,0], xy[:,1])
xyt = proj_xyz[...,:2]
ls = sgeom.LineString(xyt.tolist())
locs = axis.intersection(ls)ifnot locs:
tick =[None]else:
tick = tick_extractor(locs.xy)
_ticks.append(tick[0])# Remove ticks that aren't visible:
ticklabels = copy(ticks)whileTrue:try:
index = _ticks.index(None)except ValueError:break
_ticks.pop(index)
ticklabels.pop(index)return _ticks, ticklabels
4. 读取变量
wrfout_dir ="./data/wrfout_d02_2017-08-22_12"
wrfout = Dataset(wrfout_dir, mode="r")# 读取变量
it =8
pres = getvar(wrfout,"pressure", timeidx=it)
ua = getvar(wrfout,"ua", units="m/s", timeidx=it)
va = getvar(wrfout,"va", units="m/s", timeidx=it)
q = getvar(wrfout,"QVAPOR", timeidx=it)
times = getvar(wrfout,"times", meta=False, timeidx=ALL_TIMES)