from wrf import getvar, get_cartopy import matplotlib.pyplot as plt from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import matplotlib.transforms as mtransforms import cartopy.crs as ccrs from cartopy.io.shapereader import BasicReader import cmaps from matplotlib.path import Path from cartopy.mpl.patch import geos_to_path from mpl_toolkits.axes_grid1.inset_locator import TransformedBbox, BboxPatch, BboxConnector import shapely.geometry as sgeom from copy import copy import concurrent.futures 参考Andrew Dawson 提供的解决方法: https://gist.github.com/ajdawson/dd536f786741e987ae4e # 给出一个假定为矩形的LineString，返回对应于矩形给定边的直线。 def find_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]) # 在兰伯特投影的底部X轴上绘制刻度线 def lambert_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]) # 在兰伯特投影的左侧y轴上绘制刻度线 def lambert_yticks(ax, ticks): """Draw ricks 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]) # 获取兰伯特投影中底部X轴或左侧y轴的刻度线位置和标签 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.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) if not locs: tick = [None] else: tick = tick_extractor(locs.xy) _ticks.append(tick[0]) # Remove ticks that aren't visible: ticklabels = copy(ticks) while True: index = _ticks.index(None) except ValueError: break _ticks.pop(index) ticklabels.pop(index) return _ticks, ticklabels ## 子图连线函数 def mark_inset(parent_axes, inset_axes, loc1a=1, loc1b=1, loc2a=2, loc2b=2, **kwargs): rect = TransformedBbox(inset_axes.viewLim, parent_axes.transData) pp = BboxPatch(rect, fill=False, **kwargs) parent_axes.add_patch(pp) p1 = BboxConnector(inset_axes.bbox, rect, loc1=loc1a, loc2=loc1b, **kwargs) inset_axes.add_patch(p1) p1.set_clip_on(False) p2 = BboxConnector(inset_axes.bbox, rect, loc1=loc2a, loc2=loc2b, **kwargs) inset_axes.add_patch(p2) p2.set_clip_on(False) return pp, p1, p2 ## 等高子图设置 def add_equal_axes(ax, loc, pad, width, projection): 在原有的Axes旁新添一个等高或等宽的Axes并返回该对象. Parameters ---------- ax : Axes or array_like of Axes 原有的Axes,也可以为一组Axes构成的数组. loc : {'left', 'right', 'bottom', 'top'} 新Axes相对于旧Axes的位置. pad : float 新Axes与旧Axes的间距. width: float 当loc='left'或'right'时,width表示新Axes的宽度. 当loc='bottom'或'top'时,width表示新Axes的高度. Returns ------- ax_new : Axes 新Axes对象. # 无论ax是单个还是一组Axes,获取ax的大小位置. axes = np.atleast_1d(ax).ravel() bbox = mtransforms.Bbox.union([ax.get_position() for ax in axes]) # 决定新Axes的大小位置. if loc == 'left': x0_new = bbox.x0 - pad - width x1_new = x0_new + width y0_new = bbox.y0 y1_new = bbox.y1 elif loc == 'right': x0_new = bbox.x1 + pad x1_new = x0_new + width y0_new = bbox.y0 y1_new = bbox.y1 elif loc == 'bottom': x0_new = bbox.x0 x1_new = bbox.x1 y0_new = bbox.y0 - pad - width y1_new = y0_new + width elif loc == 'top': x0_new = bbox.x0 x1_new = bbox.x1 y0_new = bbox.y1 + pad y1_new = y0_new + width # 创建新Axes. fig = axes[0].get_figure() bbox_new = mtransforms.Bbox.from_extents(x0_new, y0_new, x1_new, y1_new) ax_new = fig.add_axes(bbox_new, projection=projection) return ax_new def process_data(): with concurrent.futures.ThreadPoolExecutor() as executor: future_d01 = executor.submit(get_data, "F:/pythonplot/geo_em.d01.nc") future_d02 = executor.submit(get_data, "F:/pythonplot/geo_em.d02.nc") # future_d03 = executor.submit(get_data, '/home/mw/input/GEO3900/geo_em.d03.nc') lon, lat, proj, hgt = future_d01.result() lon_1, lat_1,proj, hgt_1 = future_d02.result() # lon_2, lat_2,proj, hgt_2 = future_d03.result() return lon, lat, proj, hgt, lon_1, lat_1, hgt_1 def get_data(file): with Dataset(file) as f: lon = getvar(f, 'lon', meta=False) lat = getvar(f, 'lat', meta=False) proj = get_cartopy(wrfin=f) hgt = getvar(f, 'HGT_M', meta=False) hgt[hgt < 0] = 0 return lon, lat, proj, hgt lon, lat, proj, hgt, lon_1, lat_1, hgt_1 = process_data() countries = BasicReader("F:/world_adm0_Project.shp") provinces = BasicReader("F:/chinaProvince.shp") provinces1 = BasicReader("F:/gnnq.shp") cmap = cmaps.MPL_terrain #绘图部分 fig = plt.figure(figsize=(20, 12),dpi=600) plt.tight_layout() # 在figure上添加子图 ax = fig.add_axes([0.09, 0.15, 0.65, 0.8], projection=proj) #左下右上逆时针占比 contour = ax.contourf(lon, lat, hgt, transform=ccrs.PlateCarree(), cmap=cmap, levels=np.arange(0, 6500 + 500, 500)) ax.add_geometries(provinces.geometries(), linewidth=1.5, edgecolor='blue', crs=ccrs.PlateCarree(), facecolor='none') ax.add_geometries(countries.geometries(), linewidth=1.5, edgecolor='black', crs=ccrs.PlateCarree(), facecolor='none') # 添加经纬度网格和刻度 # *must* call draw in order to get the axis boundary used to add ticks: fig.canvas.draw() # Define gridline locations and draw the lines using cartopy's built-in gridliner: xticks=[80,90,100,110,120] yticks=[30,35,40,45,50] ax.gridlines(xlocs=xticks, ylocs=yticks) font3={'family':'SimHei','size':25,'color':'k'} ax.gridlines(xlocs=xticks, ylocs=yticks, draw_labels=False, linewidth=0.9, color='grey', alpha=0.6, linestyle='--') # Label the end-points of the gridlines using the custom tick makers: ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) ax.yaxis.set_major_formatter(LATITUDE_FORMATTER) lambert_xticks(ax, xticks) lambert_yticks(ax, yticks) ax.text(0.005,0.945,'D01',fontsize=24,fontweight='bold',transform=ax.transAxes) ax.tick_params(axis="x", labelsize=22) ax.tick_params(axis="y", labelsize=22) #g1.rotate_labels = False # 在d01的模拟区域上框出d02的模拟区域范围 ax.plot([lon_1[0, 0], lon_1[-1, 0]], [lat_1[0, 0], lat_1[-1, 0]], color='red', lw=2, transform=ccrs.PlateCarree()) ax.plot([lon_1[-1, 0], lon_1[-1, -1]], [lat_1[-1, 0], lat_1[-1, -1]], color='red', lw=2, transform=ccrs.PlateCarree()) ax.plot([lon_1[-1, -1], lon_1[0, -1]], [lat_1[-1, -1], lat_1[0, -1]], color='red', lw=2, transform=ccrs.PlateCarree()) ax.plot([lon_1[0, -1], lon_1[0, 0]], [lat_1[0, -1], lat_1[0, 0]], color='red', lw=2, transform=ccrs.PlateCarree()) #添加标题和颜色条 plt.title('WRF_Domain',fontsize=25, fontweight="bold", pad=0.05) cbar = plt.colorbar(contour, ax=ax, orientation='horizontal', pad=0.07, shrink=0.95) cbar.set_label('Terrain Height (m)',fontsize=25) cbar.ax.tick_params(labelsize=22) #在右上角加一个小图示例 ax_inset = add_equal_axes(ax, loc='right', pad=0.09, width=0.4, projection=proj) contour_inset = ax_inset.contourf(lon_1, lat_1, hgt_1, transform=ccrs.PlateCarree(), cmap=cmap, levels=np.arange(1000, 4500 + 500, 500)) ax_inset.add_geometries(provinces1.geometries(), linewidth=1, edgecolor='black', crs=ccrs.PlateCarree(), facecolor='none') # ax_inset.add_geometries(countries.geometries(), linewidth=1., edgecolor='black', crs=ccrs.PlateCarree(), # facecolor='none') ax_inset.text(0.005,0.945,'D02',fontsize=24,fontweight='bold',transform=ax_inset.transAxes) ax_inset.tick_params(axis="x", labelsize=22) ax_inset.tick_params(axis="y", labelsize=22) # 添加经纬度网格和刻度 fig.canvas.draw() # Define gridline locations and draw the lines using cartopy's built-in gridliner: # ax_inset.gridlines(xlocs=xticks, ylocs=yticks) xticks1=[100,101,102,103,104,105] yticks1=[37,38,39,40] ax_inset.gridlines(xlocs=xticks1, ylocs=yticks1, draw_labels=False, linewidth=0.8, color='grey', alpha=0.5, linestyle='--') # Label the end-points of the gridlines using the custom tick makers: ax_inset.xaxis.set_major_formatter(LONGITUDE_FORMATTER) ax_inset.yaxis.set_major_formatter(LATITUDE_FORMATTER) lambert_xticks(ax_inset, xticks1) lambert_yticks(ax_inset, yticks1) cbar1 = plt.colorbar(contour_inset, ax=ax_inset, orientation='vertical', pad=0.05, shrink=0.9) cbar1.ax.tick_params(labelsize=20) # cbar1.set_label('m',fontsize=25) #添加子图连线 mark_inset(ax, ax_inset, loc1a=2, loc1b=1, loc2a=3, loc2b=4, fc="none", ec='k', lw=0.85) # plt.show() plt.savefig("F:/pythonplot/wrf_insert_domain.png",bbox_inches = 'tight') #显示完整x轴图形。

参考自： WRF domain 绘制改进

原文链接： https://mp.weixin.qq.com/s/97Uu-8kfwvfVvScf3O4JGw 我通常是先根据数值天气预报的“分辨率”和要观察的经纬度范围来设计。 1、数值天气预报：      美国数值天气预报提供分辨率为1度、0.5度、1.5度的等等。     欧洲天气预报提供分辨率为0.125度、0.25度的等等。   当然是选择分辨率越小的数值天气预报越好，但这样 一直很想画垂直剖面图，借鉴了很多大佬的笔记，包括ncl和Gradas代码，最后发现还是 python 好用，然后自己根据需要改了改，此贴更多为自己记录学习下~废话不多说，附上代码。（图为T2例图，可以根据需要换参数)axe.invert_yaxis() # 翻转纵坐标。 这个环境变量，该环境变量用于设置proj4模块（Basemap模块依赖模块之一）的位置。而Anaconda中安装proj4的时候不会主动设置proj4模块的环境变量，于是就导致了上述报错。Linux系统中可将以下命令写入.bashrc文件中，source生效。 Python 脚本中，将proj4模块的位置（在conda安装目录下的。 本文主要介绍 python 对 wrf out结果文件的初步后处理操作，以及基础绘图。 wrf out后处理包括：【读取 wrf out文件、读取 wrf out文件中变量metadata及数据、对高度场进行500hPa插值、输出nc文件】 基础绘图操作包括：【设置投影和范围、绘制等值线contour和等值线标值、副高 区域 填色contourf】 仅展示初步评估 模拟 的效果，若精美绘制需要进一步的设置、细化。 如果你不喜欢命令行的操作方式，那么你可以尝试使用 python -cdo，利用 python 脚本 语言 的优势来处理气象数据。命令行的方式有其优势，比如简单易操作，可扩展性更强等，利用CDO的 python 接口也有其特有的优势，比如：通过numpy/narray可以进行直接的数据操作临时文件自动处理灵活的并行化计算条件处理操作扩展新操作符安装安装方式非常简单，运行以下命令即可：pip install cdo或... function wrf _user_getvar (nc_file, fld, it)Usage: ter = wrf _user_getvar (a, “HGT”, 0) nc_file： wrf out的nc文件名称 fld：变量名 it：时间设置。假如it=-1，则所有时间均输出。 例如：tc= wrf _user_getvar(f,"tc",-1) =========... python 中的matplotlib是一种用于创建图表的桌面绘图包（主要是2D方面）。 使用 python 对matplotlib库操作使得对图形的显现极为方便，下面是用的较多的一些用法。 建议配合I python 使用，如果通过cmd启动i python ，请使用i python –pylab启动，方便绘... wrf .getvar( wrf nc, varname, timeidx=0, method=u’cat’, squeeze=True, cache=None, meta=True, **kargs) Returns basic diagnostics from the WRF ARW model output. Below is a table of available diagnostics. Variable Name Description