WRFOUT 位温剖面和位温单格点高度图

用户11172986

发布于 2024-06-20 18:53:18

760

发布于 2024-06-20 18:53:18

前言

WRF (Weather Research and Forecasting Model) 是一种广泛用于天气预报和气候模拟的数值大气模式。通过分析WRF模型的输出数据，我们可以获得各种天气变量的空间分布及其随时间的演变情况。

在分析WRF模型输出数据时，常常需要绘制位温（Potential Temperature）剖面和位温单格点的高度图。位温是指将气块从参考高度（通常为1000 hPa）抬升或降低到某个特定高度后的温度，它是大气中的一个重要物理量，能够反映气块的垂直运动特征。

绘制位温剖面可以帮助我们理解大气的垂直结构和稳定性情况。通过观察不同高度上的位温值，我们可以推断出对流层中的温度递减率、大气边界层的稳定性等信息。而绘制位温单格点的高度图，则能够更直观地展示不同位置的位温分布及其随高度的变化趋势。

在本文中，我们将使用WRF模型的输出数据，利用Python编程语言以及相关库（如wrf-python、numpy和matplotlib）绘制位温剖面和位温单格点的高度图。我们将根据指定的站点位置和要绘制的高度层，获取相应的位温数据，并将其在图像中进行可视化展示。通过这样的绘图分析过程，我们可以更好地理解大气垂直结构以及不同位置的温度特征，为天气预报和气候研究提供有价值的参考。

位温剖面

导入库与读数据

from wrf import  to_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([os.path.join(wrfout_path, f) for f in os.listdir(wrfout_path) if f.startswith(filename_prefix)])
wrf_list = [Dataset(f) for f in wrf_files]

使用函数插值剖面

cross_start = CoordPair(lat=29, lon=105)
cross_end = CoordPair(lat=32, lon=105)

theta = getvar(wrf_list[5], 'theta')

ht = getvar(wrf_list[5], "z")

ter = getvar(wrf_list[5], "ter")
theta_cross = vertcross(theta, ht, wrfin=wrf_list[1],
                        start_point=cross_start,
                        end_point=cross_end,
                        latlon=True, meta=True)
ter_line = interpline(ter, wrfin=wrf_list[5], 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())[0])
    theta_cross_filled[0:first_idx, i] = theta_cross_filled[first_idx, i]

绘图部分(无填充版本)

# 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[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]

# 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 and  y-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[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]

# 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 and  y-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 = [30.769948959350586, 105.00334167480469]

#将经纬度坐标转换为模型坐标系（x, y）
x_y = ll_to_xy(wrf_list[5], lat_lon[0], lat_lon[1])
#读取数据
theta_s2 = getvar(wrf_list[5], "theta", timeidx=0)[:, x_y[0], x_y[1]] * units.degC
h = getvar(wrf_list[5], "height", timeidx=0)[:, x_y[0], x_y[1]] * 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)