WRFOUT 位温剖面和位温单格点高度图
前言
WRF (Weather Research and Forecasting Model) 是一种广泛用于天气预报和气候模拟的数值大气模式。通过分析WRF模型的输出数据，我们可以获得各种天气变量的空间分布及其随时间的演变情况。
在分析WRF模型输出数据时，常常需要绘制位温（Potential Temperature）剖面和位温单格点的高度图。位温是指将气块从参考高度（通常为1000 hPa）抬升或降低到某个特定高度后的温度，它是大气中的一个重要物理量，能够反映气块的垂直运动特征。
绘制位温剖面可以帮助我们理解大气的垂直结构和稳定性情况。通过观察不同高度上的位温值，我们可以推断出对流层中的温度递减率、大气边界层的稳定性等信息。而绘制位温单格点的高度图，则能够更直观地展示不同位置的位温分布及其随高度的变化趋势。
在本文中，我们将使用WRF模型的输出数据，利用Python编程语言以及相关库（如wrf-python、numpy和matplotlib）绘制位温剖面和位温单格点的高度图。我们将根据指定的站点位置和要绘制的高度层，获取相应的位温数据，并将其在图像中进行可视化展示。通过这样的绘图分析过程，我们可以更好地理解大气垂直结构以及不同位置的温度特征，为天气预报和气候研究提供有价值的参考。
位温剖面
导入库与读数据
from wrf importto_np, getvar, interplevel, get_cartopy, latlon_coords,ALL_TIMES,CoordPair,vertcross,interpline,ll_to_xy
import numpy as np
from netCDF4 import Dataset
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cmaps
from glob import glob
import metpy.calc as mpcalc
import metpy.constants as constants
import os
# 定义 WRF 文件夹路径和文件名前缀
wrfout_path = "/home/mw/input/wrfout3385/"
filename_prefix = "wrfout_d02_"
wrf_files = sorted()
wrf_list = 使用函数插值剖面
cross_start = CoordPair(lat=29, lon=105)
cross_end = CoordPair(lat=32, lon=105)
theta = getvar(wrf_list, 'theta')
ht = getvar(wrf_list, "z")
ter = getvar(wrf_list, "ter")
theta_cross = vertcross(theta, ht, wrfin=wrf_list,
start_point=cross_start,
end_point=cross_end,
latlon=True, meta=True)
ter_line = interpline(ter, wrfin=wrf_list, start_point=cross_start,
end_point=cross_end)
# Get the height coordinate values
xs = np.arange(0, theta_cross.shape[-1], 1)
ys = to_np(theta_cross.coords["vertical"])
填充
theta_cross_filled = np.ma.copy(to_np(theta_cross))
for i in range(theta_cross_filled.shape[-1]):
column_vals = theta_cross_filled[:,i]
first_idx = int(np.transpose((column_vals > -200).nonzero()))
theta_cross_filled = theta_cross_filled
绘图部分(无填充版本)
# Get the height coordinate values
xs = np.arange(0, theta_cross.shape[-1], 1)
ys = to_np(theta_cross.coords["vertical"])
# Create the figure
fig, ax = plt.subplots(figsize=(8, 6))
# Plot the theta cross section
contours = ax.contourf(xs,ys,
to_np(theta_cross),
levels=np.arange(300,500,10),
cmap=cmaps.temp_19lev)
# Fill in the mountain area
ht_fill = ax.fill_between(xs, 0, to_np(ter_line),
facecolor="saddlebrown")
# Set the x-ticks to use latitude and longitude labels
coord_pairs = to_np(theta_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape)
x_labels =
# Set the desired number of x ticks below
num_ticks = 5
thin = int((len(x_ticks) / num_ticks) + .5)
ax.set_xticks(x_ticks[::thin])
ax.set_xticklabels(x_labels[::thin], rotation=45, fontsize=8)
# Set the x-axis andy-axis labels
ax.set_xlabel("Latitude, Longitude", fontsize=12)
ax.set_ylabel("Height (m)", fontsize=12)
# Add a title
ax.set_title("Cross-Section of Potential Temperature", fontsize=14)
plt.show()
/opt/conda/lib/python3.7/site-packages/cmaps/cmaps.py:3869: UserWarning: Trying to register the cmap 'temp_19lev' which already exists.
matplotlib.cm.register_cmap(name=cname, cmap=cmap)
绘图部分（填充版本）
# Get the height coordinate values
xs = np.arange(0, theta_cross.shape[-1], 1)
ys = to_np(theta_cross.coords["vertical"])
# Create the figure
fig, ax = plt.subplots(figsize=(8, 6))
# Plot the theta cross section
contours = ax.contourf(xs,ys,
to_np(theta_cross_filled),
levels=np.arange(300,500,10),
cmap=cmaps.temp_19lev)
# Fill in the mountain area
ht_fill = ax.fill_between(xs, 0, to_np(ter_line),
facecolor="saddlebrown")
# Set the x-ticks to use latitude and longitude labels
coord_pairs = to_np(theta_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape)
x_labels =
# Set the desired number of x ticks below
num_ticks = 5
thin = int((len(x_ticks) / num_ticks) + .5)
ax.set_xticks(x_ticks[::thin])
ax.set_xticklabels(x_labels[::thin], rotation=45, fontsize=8)
# Set the x-axis andy-axis labels
ax.set_xlabel("Latitude, Longitude", fontsize=12)
ax.set_ylabel("Height (m)", fontsize=12)
# Add a title
ax.set_title("Cross-Section of Potential Temperature", fontsize=14)
plt.show()
位温格点高度图
方法一
theta_s = theta_cross[:,20]
theta_s
可以看出从剖面取出的值还是有对应的经纬度，可直接绘制单格点高度分布
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(theta_s,ys)
ax.set_xlabel("theta(K) at lon 105.00 lat 30.77 ", fontsize=12)
ax.set_ylabel("Height (m)", fontsize=12)
Text(0, 0.5, 'Height (m)')
方法二
#指定要提取的经纬度坐标点
lat_lon =
#将经纬度坐标转换为模型坐标系（x, y）
x_y = ll_to_xy(wrf_list, lat_lon, lat_lon)
#读取数据
theta_s2 = getvar(wrf_list, "theta", timeidx=0)[:, x_y, x_y] * units.degC
h = getvar(wrf_list, "height", timeidx=0)[:, x_y, x_y] * units.degC
# plot
fig, ax2 = plt.subplots(figsize=(8, 6))
ax2.plot(theta_s2,h)
ax2.set_xlabel("theta(K) at lon 105.00 lat 30.77 ", fontsize=12)
ax2.set_ylabel("Height (m)", fontsize=12)
Text(0, 0.5, 'Height (m)')
差别
乍一看没差别，实际上两者的shape是不同的theta_s.shape,theta_s2.shape
((100,), (49,))
剖面取的点是插值之后的，层数达100，而直接取的单格点仅有模式设置的49层。
这单点高度图只是随手之作，大家有更好的办法可以评论区讨论讨论。
小结
做得仓促，没有细化绘图。从剖面再取格点貌似绕了远路（难道我会告诉你只是剖面图的副产物吗）
这时候有同学要问了，这地形图怎么这么难看啊？都说是仓促作图。废话少说赶紧点赞。
完整代码与文件在此
文章来源于微信公众号：气ython风雨
页:
[1]