首页 > 解决方案 > 如何将等压相对湿度数据转换为二维数组形式?

问题描述

我正在尝试在类似于MetPy xarray 教程中的恒定压力水平(500hPa)下创建相对湿度的等高线图。我已经使用 Siphon 包获取了数据,并将其解析为一个似乎是二维的数组,时间和高度固定,纬度/经度变化:

<xarray.DataArray 'relative_humidity' (time: 1, lat: 141, lon: 121)>
array([[[100. , 100. , ...,  48.5,  48.1],
        [100. , 100. , ...,  42.8,  41.1],
        ...,
        [  9.5,   9.4, ...,  20.7,  18.7],
        [  9.5,   9.9, ...,  23.8,  21.1]]], dtype=float32)
Coordinates:
    reftime   (time) datetime64[ns] 2020-03-24T12:00:00
  * time      (time) datetime64[ns] 2020-03-24T18:00:00
    isobaric  float32 50000.0
  * lat       (lat) float32 55.0 54.75 54.5 54.25 54.0 ... 20.75 20.5 20.25 20.0
  * lon       (lon) float32 270.0 270.25 270.5 270.75 ... 299.5 299.75 300.0
    crs       object Projection: latitude_longitude

为了得到这个数组,我使用了代码:

data = ncss.get_data(query)

#Parse data using MetPy
ds = xr.open_dataset(NetCDF4DataStore(data))
data = ds.metpy.parse_cf()

#Rename variables to useful things
data = data.rename({
    'Vertical_velocity_pressure_isobaric': 'omega',
    'Relative_humidity_isobaric': 'relative_humidity',
    'Temperature_isobaric': 'temperature',
    'u-component_of_wind_isobaric': 'u',
    'v-component_of_wind_isobaric': 'v',
    'Geopotential_height_isobaric': 'height'
})

#Get data specific to 500mb
zH5 = data['height'].metpy.sel(vertical=850 * units.hPa)
zH5_crs = zH5.metpy.cartopy_crs

#Define coordinates
vertical, = data['temperature'].metpy.coordinates('vertical')
time = data['temperature'].metpy.time
x, y = data['height'].metpy.coordinates('x', 'y')
lat, lon = xr.broadcast(y, x)

#Create relative humidity array
rel_hum = data['relative_humidity'].metpy.sel(vertical=500*units.hPa)

但是,当我将 RH 绘制为填充轮廓时:

rh = ax.contourf(x, y, rel_hum, levels=[70, 80, 90, 100], colors=['#99ff00', '#00ff00', '#00cc00'])

我收到一条错误消息:

TypeError: Input z must be 2D, not 3D

在阅读了关于类似问题的另一篇 SO 帖子后,我明白为什么轮廓函数的第三个参数需要是一个二维数组,但我不确定为什么我的程序在这里(它非常模仿来自MetPy docs ) 不会产生这样一个能够被绘制的数组。

标签: python-xarraymetpy

解决方案


因此,虽然您的数据似乎是 2D 的,因为您只有一次,但基于此输出:

<xarray.DataArray 'relative_humidity' (time: 1, lat: 141, lon: 121)>

数据仍被视为 3D,因为您仍然有时间维度(即使它的大小为 1)。一种方法是修改您的呼叫以.sel选择特定时间(或添加另一个呼叫以.sel执行此操作)。

我首选的解决方法是调用.squeeze(). squeeze()是一种删除所有 size-1 维度的方法,解决了这个非常常见的问题,即具有一些无关维度(实际上并没有添加到数组中的元素总数中):

rel_hum = data['relative_humidity'].metpy.sel(vertical=500*units.hPa).squeeze()

这也有一个效果,如果你真的有 3D 数据,它不会改变形状。根据您的应用程序,这可能是也可能不是一件好事。(如果不好,使用.sel选择时间是更好的选择。)


推荐阅读