你怎么知道我深蹲100斤了？学习笔记 | metpy插值气压层时需要注意什么

用户11172986

发布于 2026-04-24 18:54:56

780

metpy 插值到等压面总是报错？90 %的人都忽略了单位这件事——一次 Xarray 维度踩坑全记录

环境设置

安装依赖

!pip install metpy

正文

近来又潮又热，去江苏出差的朋友感慨江苏怎么还没北京潮，还有读者在评论区发了北京已有广东特有的美洲大蠊，真是奇哉怪哉。

言归正传，如果你在气象数据处理圈子里混，多多少少肯定碰过这段代码。没碰过也没关系，现在碰过了—— 把模式层数据插值到指定等压面：

import metpy.interpolate as mpint  

plevs = np.array([1000, 850, 700, 500])*100    # Pa  
new_u = mpint.log_interpolate_1d(  
            plevs, pressure, u, axis=1)

很多时候，跑得好好的脚本，某天突然蹦出一段让人摸不着头脑的报错：

TypeError: no implementation found for 'numpy.apply_along_axis'

或者：

ValueError: Cannot apply units on ...

别慌，这其实不是 bug，而是 90 % 的人都没有掌握的单位+维度“潜规则”。实际上就是metpy老掉牙的单位问题。今天咱们一次性说清楚。

01 场景复现：到底哪一步炸裂？

官方 issue #1889 给了一个最简单、却最容易忽略的复现：

Xarray DataArray ➜ 4 维（time, lev, lat, lon）
Numpy ndarray ➜ 同样形状

前者直接炸；后者顺溜过关。

这是为什么呢

02 背后真相：Xarray 带来的隐藏属性

输入类型	has units?	has coords?	apply_along_axis 是否能识别
ndarray	✗	✗	✅ 直接上裸数组
DataArray	maybe	✅	⚠️ metpy 必须拆 unit，否则 numpy 拒接

metpy 为了保持单位一致性，会在内部把 pressure 先变成 pint.Quantity。但 new_levels 是裸 numpy，两者单位不统一 → _strip_matching_units() 拆不掉 → 最后一步 log_interpolate_1d 掉到 apply_along_axis，numpy 见到带 unit 的数组直接罢工。

03 三分钟排查清单

所有输入都要有单位

from metpy.units import units  
plevs = (np.array([...,...,...])*100) * units.Pa   # 千万别省 * units

如果懒得写单位？那就全部无量纲

plevs = plevs * units('')  
pressure = pressure.metpy.quantify()

用 .metpy.dequantify() 随时把数据剥干净仅当你确定后续不再依赖单位，或者喂给其它裸 numpy 库。

维度顺序要一致

pressure = pressure.transpose(*u.dims)

04 实战示例

import xarray as xr
import numpy as np
from metpy.interpolate import log_interpolate_1d
from metpy.units import units

# 读取数据并确保单位系统
ds = xr.open_dataset('atmos.nc', decode_times=False)
print(ds)

<xarray.Dataset> Size: 58MB
Dimensions:      (ilev: 31, lev: 30, lat: 64, lon: 128, time: 1)
Coordinates:
  * lev          (lev) float32 120B 3.643 7.595 14.36 ... 957.5 976.3 992.6
  * ilev         (ilev) float32 124B 2.255 5.032 10.16 ... 967.5 985.1 1e+03
  * lat          (lat) float32 256B -87.86 -85.1 -82.31 ... 82.31 85.1 87.86
  * lon          (lon) float32 512B 0.0 2.812 5.625 8.438 ... 351.6 354.4 357.2
  * time         (time) float64 8B 791.0
Data variables: (12/130)
    P0           float32 4B ...
    ntrm         int32 4B ...
    ntrn         int32 4B ...
    ntrk         int32 4B ...
    ndbase       int32 4B ...
    nsbase       int32 4B ...
    ...           ...
    TAUX         (time, lat, lon) float32 33kB ...
    TAUY         (time, lat, lon) float32 33kB ...
    SNOWH        (time, lat, lon) float32 33kB ...
    UTGW         (time, lev, lat, lon) float32 983kB ...
    VTGW         (time, lev, lat, lon) float32 983kB ...
    VVPUU        (time, lev, lat, lon) float32 983kB ...
Attributes:
    Conventions:            NCAR-CSM
    source:                 Data converted from CCM History Tape Format
    case:                   prgcld01
    title:                  Phils prognostic cld water + tendency physics code
    hybrid_sigma_pressure:  \nPressure at a grid point (lon(i),lat(j),lev(k))...
    history:                \nMon May 17 12:31:49 1999>


u = ds.U * units('m/s')  # 强制添加单位（如果原始数据没有）
ps = ds.PS * units.Pa    # 表面气压必须带单位

# 系数处理（CAM模式常用）
hyam = ds.hyam * units('dimensionless')
hybm = ds.hybm * units('dimensionless')
p0 = 1000 * units.hPa

# 计算全层气压
pressure = hyam * p0 + hybm * ps
pressure = pressure.transpose(*u.dims) 

# 创建等压面层级（hPa转Pa）
new_levels = np.array([1000, 950, 800, 700, 600, 500, 400,200]) * units.hPa
# 对数插值
u_isobaric = log_interpolate_1d(new_levels, pressure, u.metpy.quantify(), axis=1)
# 验证结果
print("插值结果单位:", u_isobaric.units)  # 应该显示 "meter / second"
print("NaN值数量:", np.isnan(u_isobaric).sum())  # 检查异常值

插值结果单位: meter / second
NaN值数量: 7466

print(u_isobaric[0,0])

[[nan nan nan ... nan nan nan] [nan nan nan ... nan nan nan] [nan nan nan ... nan nan nan] ... [-2.44014092908361 -2.5727052585667027 -2.6622389092978307 ... nan  -2.010691574429729 -2.256805604584203] [-2.30256332842527 -2.2468413152222126 -2.181935846149208 ...  -2.3940958331870448 -2.3783463214362417 -2.3471175183966806] [-2.0253718473951325 -1.9608527770423514 -1.8944898069814329 ...  -2.203676466735914 -2.1471989814030255 -2.0876272126679476]] meter / second