请教：怎么把ERA5中臭氧质量混合比转换成质量浓度

diva211

想画一个臭氧的垂直分布-时间剖面图，并且单位以微克每立方米表示，在单位转换的步骤会报错
f = addfile("k:/MI/ozone.nc")
udata = f['u'][:,::-1,"36.66","117.04"]
vdata = f['v'][:,::-1,"36.66","117.04"]
rmdata=f['r'][:,::-1,"36.66","117.04"]
tmdata=f['t'][:,::-1,"36.66","117.04"]
Tmdata=f['valid_time'][:]
wmdata=f['w'][:,::-1,"36.65","117.04"]
omdata = f['o3'][:, ::-1, "36.65", "117.15"]
u=udata.T*2.5
v=vdata.T*2.5
r=rmdata.T
t=tmdata.T
w=wmdata.T
o=omdata.T


# 获取气压层数据
level = o.dimvalue(0)
level = level[::-1]  # 反转气压层顺序


# 将气压层转换为标准高度
height = meteolib.pressure_to_height_std(level)
height = height[:] / 1000  # 转换为千米
o.setdimvalue(0, height)


# 获取时间维度数据
tt = o.dimvalue(1)
#世界时转北京时
tt = tt + 3.6 * 8

# 计算空气密度 (kg/m³)
Rd = 287.05  # 干空气气体常数 (J/(kg·K))
air_density = (level * 100) / (Rd * t)  # 气压转换为Pa，计算密度


# 将臭氧质量混合比转换为质量浓度 (µg/m³)
ozone_concentration = o * air_density * 1e9  # 1e9 将 kg/m³ 转换为 µg/m³


# 绘制臭氧质量浓度
layer1 = contourf(tt, height, ozone_concentration)
colorbar(layer1, fontsize=16, bold="True")


# 设置x轴
xaxis(axistype='time', timetickformat='ddHH')
xticks(fontsize=15, bold=True)
xlabel('Time (UTC)', fontsize=15, bold="True")
xlim(tt[169], tt[288])


# 设置y轴
ylim(height.min(), 1.458)
levy = array([1000, 975, 950, 925, 900, 875, 850])
yticks(o.dimvalue(0), levy)
ylabel('Pressure (hPa)', fontsize=15, bold="True")

MeteoInfo

哪行代码出错？具体的错误信息是什么？

diva211

MeteoInfo 发表于 2025-2-18 11:27
哪行代码出错？具体的错误信息是什么？

>>> run script...
Traceback (most recent call last):
  File "<iostream>", line 37, in <module>
  File "K:\MeteoInfo\pylib\mipylib\numeric\core\dimarray.py", line 341, in __rdiv__
    r = super(DimArray, self).__rdiv__(other)
  File "K:\MeteoInfo\pylib\mipylib\numeric\core\_ndarray.py", line 391, in __rdiv__
    raise ValueError('Dimension missmatch, can not broadcast!')
ValueError: Dimension missmatch, can not broadcast!

王老师，这是报错信息，报错行已标红

MeteoInfo

diva211 发表于 2025-2-18 14:10
>>> run script...
Traceback (most recent call last):
  File "", line 37, in

你试试先不要进行数组转置，计算结束后画图前再转置

diva211

MeteoInfo 发表于 2025-2-18 14:55
你试试先不要进行数组转置，计算结束后画图前再转置

按照您说的先计算再转置，成功了，非常感谢！

f = addfile("k:/MI/ozone.nc")
udata = f['u'][:,::-1,"36.66","117.04"]
vdata = f['v'][:,::-1,"36.66","117.04"]
rdata=f['r'][:,::-1,"36.66","117.04"]
tdata=f['t'][:,:,"36.66","117.04"]
Tdata=f['valid_time'][:]
wdata=f['w'][:,::-1,"36.65","117.04"]
odata = f['o3'][:,:, "36.65", "117.15"]

# 获取气压层数据
level = odata.dimvalue(1)
# 计算空气密度 (kg/m&#179;)
Rd = 287.05  # 干空气气体常数 (J/(kg·K))
air_density = (level * 100) / (Rd * tdata)  # 气压转换为Pa，计算密度
# 将臭氧质量混合比转换为质量浓度 (&#181;g/m&#179;)
ozone_concentration = odata * air_density * 1e9  # 1e9 将 kg/m&#179; 转换为 &#181;g/m&#179;
oc=ozone_concentration[:,::-1]
o3=oc.T

tdata=f['t'][:,::-1,"36.66","117.04"]
odata = f['o3'][:,::-1, "36.65", "117.15"]
u=udata.T*2.5
v=vdata.T*2.5
r=rdata.T
t=tdata.T
w=wdata.T
o=odata.T


# 将气压层转换为标准高度
level = o.dimvalue(0)
level = level[::-1]  # 反转气压层顺序
height = meteolib.pressure_to_height_std(level)
height = height[:] / 1000  # 转换为千米
o.setdimvalue(0, height)

# 获取时间维度数据
tt = o.dimvalue(1)
#世界时转北京时
tt = tt + 3.6 * 8

# 绘制臭氧质量浓度
layer1 = contourf(tt, height, o3)
colorbar(layer1, fontsize=16, bold="True")

# 设置x轴
xaxis(axistype='time', timetickformat='ddHH')
xticks(fontsize=15, bold=True)
xlabel('Time (UTC)', fontsize=15, bold="True")
xlim(tt[169], tt[288])

# 设置y轴
ylim(height.min(), 1.6)
levy = array([1000, 975, 950, 925, 900, 875, 850])
yticks(o.dimvalue(0), levy)
ylabel('Pressure (hPa)', fontsize=15, bold="True")