当前位置：首页>python>Python | MCS追踪 |tobac

Python | MCS追踪 |tobac

2026-09-11 02:25:16

MCS-tracking based on ICON healpix grid

准备数据

降水、OLR、表面气温

cat = intake.open_catalog('https://data.nextgems-h2020.eu/catalog.yaml')

pd.DataFrame(cat.ICON.C5.AMIP_CNTL.describe()["user_parameters"])
icon_ds = cat.ICON.C5.AMIP_CNTL(zoom=8, time="PT20M", chunks={}).to_dask().pipe(egh.attach_coords)
icon_ds
icon_ds = icon_ds[["pr", "rlut", "tas"]].sel(time=slice(datetime(1980,1,1), datetime(1980,1,2))).load()
img = egh.healpix_show(icon_ds.rlut[42])
plt.colorbar(img, orientation="horizontal", label="ToA OLR [W m-2]")
plt.title("OLR")

OLR常用于追踪砧状云，但在全球模拟中，OLR随纬度的变化会导致冬季半球极地附近出现误报。降水也可以用于追踪强降水风暴，但会错过砧状云的衰减阶段。

为了解决 OLR 的问题，我们可以用地表气温的黑体辐射减去 OLR 来近似大气对 OLR 的贡献。这样得到的场中，高云的数值较大，而寒冷晴朗的天空场景的数值较小。

surface_olr_diff = (icon_ds.tas**4 * 5.67e-8) - icon_ds.rlut

可以使用以下函数根据 OLR 近似计算 11 微米亮度温度 (BT)：

defget_tb(olr):
"""
    This function converts outgoing longwave radiation to brightness temperatures.

    Args:
        olr(xr.DataArray or numpy array): 2D field of model output with OLR

    Returns:
        tb(xr.DataArray or numpy array): 2D field with estimated brightness temperatures
    """
# constants
    aa = 1.228
    bb = -1.106e-3# K−1
# Planck constant
    sigma = 5.670374419e-8# W⋅m−2⋅K−4

# flux equivalent brightness temperature
    Tf = (abs(olr) / sigma) ** (1.0 / 4)
    tb = (((aa**2 + 4 * bb * Tf) ** (1.0 / 2)) - aa) / (2 * bb)
return tb

将 ICON 数据转换为常规经纬度网格以进行跟踪

最简单的网格重划分方法是，在找到与经纬度网格匹配的单元格后，使用最近邻healpix.ang2pix算法.

deflatlon_remap_cells(grid_spacing, nside, order="nest"):
    lon = xr.DataArray(
        np.arange(grid_spacing / 2, 360, grid_spacing),
        dims=("lon",),
        name="lon",
        attrs=dict(units="degrees", standard_name="longitude")
    )

    lat = xr.DataArray(
        np.arange(90 - grid_spacing / 2, -90, -grid_spacing),
        dims=("lat",),
        name="lat",
        attrs=dict(units="degrees", standard_name="latitude")
    )

    pix = healpix.ang2pix(
        nside,
        *np.meshgrid(lon.values, lat.values),
        nest=(order == "nest"),
        lonlat=True
    )

return xr.DataArray(
        pix,
        coords={"lat": lat, "lon": lon},
        dims=("lat", "lon"),
        name="cell"
    )

使用 xarray 的最近邻选择功能进行重新网格化：

zoom = 8
nside = 2 ** zoom
grid_spacing = 64 / 2**zoom
dxy = 1e5 * grid_spacing
grid_spacing
nside
pix = latlon_remap_cells(grid_spacing, nside, order="nest")

选择变量进行网格转换：

surface_olr_diff = (icon_ds.tas**4 * 5.67e-8) - icon_ds.rlut
surface_olr_diff_gridded = surface_olr_diff.drop_vars(["lat", "lon"], errors="ignore").sel(cell=pix, method="nearest")
surface_olr_diff_gridded

进行特征检测和跟踪

_目标检测过程在tobac_中包括四个步骤：

特征检测
追踪
合并/拆分
分割

首先，我们

使用多个阈值和 x 方向上的周期性边界进行特征检测。

然后我们使用运动预测进行跟踪：

合并/分割操作会找到起始/结束距离小于 dxy*4 且与其他特征相距 1 个时间步的特征：

Features = tobac.feature_detection_multithreshold(

surface_olr_diff_gridded.drop_vars(["cell", "crs"], errors="ignore"),

dxy=dxy, threshold=[250, 275, 300, 325], min_distance=dxy*4, PBC_flag="hdim_2"

)

Tracks = tobac.linking_trackpy(

Features, None, dt=900, dxy=dxy, v_max=50,

method_linking="predict", stubs=2,

PBC_flag="hdim_2", min_h2=0, max_h2=surface_olr_diff_gridded.lon.size

)

merges_splits = tobac.merge_split.merge_split_MEST(

Tracks, dxy=dxy, distance=dxy*4, frame_len=1, PBC_flag="hdim_2", min_h2=0, max_h2=surface_olr_diff_gridded.lon.size

)

merges_splits

最后，我们进行分割，以找到每个特征周围阈值为 225 W 的区域。

segmentation_mask, Tracks = tobac.segmentation.segmentation(

Tracks, surface_olr_diff_gridded, dxy=dxy, threshold=225, PBC_flag="hdim_2", seed_3D_flag="box"

)

Tracks