如何从4维dataarray中提取对角等非水平方向的特定区域

要提取任意非水平方向的剖面/区域，可以通过坐标插值+自定义坐标序列的方式实现，具体操作如下：

步骤1：定义目标区域的坐标序列 #

首先确定你要提取的区域的坐标点，以对角斜线为例，先定义两个端点的x、y坐标，再生成连续的插值点：

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

# 自定义对角剖面的两个端点坐标，可根据你的需求修改为任意位置
x_start, y_start = -200000, -200000
x_end, y_end = 150000, 200000

# 生成两点间100个均匀分布的坐标点，点数越多剖面越平滑
n_points = 100
x_profile = np.linspace(x_start, x_end, n_points)
y_profile = np.linspace(y_start, y_end, n_points)

# 可选：计算每个点到起点的距离，用于后续可视化的横轴
dist_km = np.sqrt((x_profile - x_start)**2 + (y_profile - y_start)**2) / 1000

步骤2：提取对应坐标的数值 #

用xarray自带的插值/索引方法提取目标区域的数值：

提取固定高度层的斜向剖面（1维序列）

对应你需求里time=3、z=4的固定层：

# 线性插值法提取，更平滑
ref_profile = tkt16l['REF'].interp(
    time=3, 
    z=4, 
    x=xr.DataArray(x_profile, dims='profile'), 
    y=xr.DataArray(y_profile, dims='profile')
)

# 若不需要插值，直接取最邻近的网格点，用以下代码即可
# ref_profile = tkt16l['REF'].sel(
#     time=3, 
#     z=4, 
#     x=xr.DataArray(x_profile, dims='profile'), 
#     y=xr.DataArray(y_profile, dims='profile'),
#     method='nearest'
# )

提取沿斜向的垂直剖面（2维，包含z方向所有层）

如果需要看这条斜线上不同高度的反射率分布，不用固定z维度即可：

ref_vert_profile = tkt16l['REF'].interp(
    time=3, 
    x=xr.DataArray(x_profile, dims='profile'), 
    y=xr.DataArray(y_profile, dims='profile')
)

步骤3：可视化结果 #

# 可视化固定高度的斜向剖面
plt.figure(figsize=[10,3], dpi=100)
plt.plot(dist_km, ref_profile, c='k')
plt.xlabel('距离起点的距离(km)')
plt.ylabel('反射率(dBZ)')
plt.grid(alpha=0.3)
plt.show()

# 可视化沿斜线的垂直剖面
plt.figure(figsize=[10,4], dpi=100)
ref_vert_profile.plot.contourf(cmap='gray_r', levels=np.linspace(25,60,5))
plt.xticks(ticks=np.arange(0, n_points, 20), labels=np.round(dist_km[::20], 1))
plt.xlabel('距离起点的距离(km)')
plt.ylabel('高度(m)')
plt.show()

# 如果要在原有水平切片上标记你提取的剖面位置，在你原来的绘图代码后加一行即可
# tkt16l['REF'][3,4].plot.contourf(cmap='gray_r', levels=np.linspace(25,60,5), xlim=[-200000,150000],ylim=[-200000,200000])
# plt.plot(x_profile, y_profile, c='red', linewidth=2, label='提取的剖面')
# plt.legend()

注意事项 #

• 上述代码默认你的DataArray维度名和坐标名为time/z/y/x，如果你的数据集命名不同，对应修改即可。
• 如果要提取折线、多边形等其他形状的区域，只要把x_profile和y_profile替换为对应形状的坐标点序列即可，多边形区域可以结合xarray.where做掩膜提取。

内容的提问来源于stack exchange，提问作者Wind Circulation