MCS-tracking based on ICON healpix grid
准备数据
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
image进行特征检测和跟踪
_目标检测过程在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
可以将其网格再转化回去，进行绘图分析，但是我更熟悉lat-lon的网格,先记录到这吧，后面绘图分析的跑完再看。
★https://tobac.readthedocs.io/en/latest/examples/Example_ICON_MCS_tracking/Example_ICON_MCS_tracking.html
