代码实战 | 基于era5的温度平流计算

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy.units import units
import metpy.calc as mpcalc

def plot_temperature_advection(lon, lat, data, title):
    fig = plt.figure(figsize=(12, 6))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    ax.add_feature(cfeature.COASTLINE)
    cs = ax.contourf(lon, lat, data, levels=10, cmap='RdBu_r', transform=ccrs.PlateCarree())
    plt.colorbar(cs, orientation='vertical', label='Temperature Advection (K/s)')
    ax.set_title(title)
    plt.show()

def metpy_method(lon, lat, u, v, t):
    # 转换为MetPy单位数组
    t = t * units.K
    u = u * units('m/s')
    v = v * units('m/s')
    lon = lon * units.degrees
    lat = lat * units.degrees
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
    # 计算温度平流
    adv = mpcalc.advection(t, u=u, v=v, dx=dx, dy=dy)
    return adv.magnitude  # 返回无单位的numpy数组

# 读取ERA5数据
def load_era5_data(filepath):
    # 打开netCDF文件
    ds = xr.open_dataset(filepath)
    ds = ds.isel(time=0, level=25)  # 示例：第一个时间步和第一个层面
    print(ds)
    u = ds['u'].values  # 纬向风 (m/s)
    v = ds['v'].values  # 经向风 (m/s)
    t = ds['t'].values  # 温度 (K)
    lon = ds['longitude'].values
    lat = ds['latitude'].values
    
    return lon, lat, u, v, t

# 使用实际ERA5数据路径
era5_file = '/home/mw/input/era58091/ERA5-2023-08_pl.nc'
lon0, lat0, u0, v0, t = load_era5_data(era5_file)
temadvection_metpy = metpy_method(lon0, lat0, u0, v0, t)
plot_temperature_advection(lon0, lat0, temadvection_metpy, "Temperature Advection (MetPy)")

<xarray.Dataset> Size: 2MB
Dimensions:    (longitude: 361, latitude: 241)
Coordinates:
  * longitude  (longitude) float32 1kB 90.0 90.25 90.5 ... 179.5 179.8 180.0
  * latitude   (latitude) float32 964B 70.0 69.75 69.5 69.25 ... 10.5 10.25 10.0
    level      int32 4B 700
    time       datetime64[ns] 8B 2023-08-02
Data variables:
    z          (latitude, longitude) float32 348kB ...
    r          (latitude, longitude) float32 348kB ...
    t          (latitude, longitude) float32 348kB ...
    q          (latitude, longitude) float32 348kB ...
    u          (latitude, longitude) float32 348kB ...
    v          (latitude, longitude) float32 348kB ...
    w          (latitude, longitude) float32 348kB ...
Attributes:
    Conventions:  CF-1.6
    history:      2024-01-09 16:38:26 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

def correct_method(lon, lat, u, v, t):
    a = 6371000.0# 地球半径，单位米
    deg2rad = np.pi / 180
    
    # 转换为numpy数组
    lon = np.asarray(lon)
    lat = np.asarray(lat)
    u = np.asarray(u)
    v = np.asarray(v)
    t = np.asarray(t)
    
    # 计算网格间距（考虑地球曲率）
    dlon = np.gradient(lon) * deg2rad  # 转为弧度
    dlat = np.gradient(lat) * deg2rad  # 转为弧度
    
    # 转换为二维数组便于广播
    dlon = dlon.reshape((1, -1))
    dlat = dlat.reshape((-1, 1))
    coslat = np.cos(lat * deg2rad).reshape((-1, 1))
    # 计算实际距离（米）
    dx = a * coslat * dlon
    dy = a * dlat
    # 计算温度梯度（使用中心差分）
    dT_dx = np.gradient(t, axis=1) / dx
    dT_dy = np.gradient(t, axis=0) / dy
    # 温度平流（-V·∇T）
    temadvection = -(u * dT_dx + v * dT_dy)
    
    return temadvection

temadvection_correct = correct_method(lon0, lat0, u0, v0, t)
plot_temperature_advection(lon0, lat0, temadvection_correct, "Correct Temperature Advection Calculation")

# 比较方法的结果
print("正确方法结果范围:", np.nanmin(temadvection_correct ), np.nanmax(temadvection_correct ))
print("MetPy方法结果范围:", np.nanmin(temadvection_metpy), np.nanmax(temadvection_metpy))

# 计算MetPy与正确方法的差异
diff = temadvection_metpy - temadvection_correct 
print("MetPy与正确方法的最大差异:", np.nanmax(np.abs(diff)))

正确方法结果范围: -0.000914788 0.0010622571
MetPy方法结果范围: -0.0009147886257277709 0.001062257697235703
MetPy与正确方法的最大差异: 8.654054701647786e-05