警告

本文仅分享处理思路，供参考和学习，只有部分代码片段。请遵守中国气象数据网的使用条款，合理使用数据。

上一篇文章确定了中国气象数据网•综合实况的 RADAR_L3_MST_CREF 雷达拼图 BIN 瓦片是 256x256、每像素 2 字节的大端整数，但坐标顺序处理不当，绘制结果细节有问题。这篇文章记录后续的处理：确定瓦片网格公式、验证坐标方向、处理缺档瓦片，并将结果拼接为带经纬度坐标的 NetCDF 数据集。

瓦片网格

瓦片 URL 形如：

1

https://image.data.cma.cn/tiles/China/{product}/{YYYYMMDD}/{HH}/{MM}/{z}/{y}/{x}.bin

路径中的时间为北京时（UTC+8），6 分钟一档（分钟取 00/06/…/54），发布滞后约 20–30 分钟。

该网站前端加载的 canvasLeflet chunk 里包含 tileActualBounds 的计算逻辑，解混淆后得到的是等经纬度（plate carrée）网格：

1
2
3

degPerTile = 360 / 2^z
西边界 lon = x * degPerTile - 180
南边界 lat = y * degPerTile - 90     # y 自南向北计数

纬度方向每块瓦片同样跨 360/2^z 度，因此只有 y < 2^(z-1) 的行有意义。可用层级为 z=2 到 z=5，z=5 时单像素约 360/(256*32) ≈ 0.044°。

若按标准 Web Mercator XYZ 解读，经度方向的公式恰好一致（列公式相同），纬度方向则会整体错位、行序错乱。验证方法是叠加国境线后检查孤立的雷达覆盖圆盘：单个雷达站的覆盖范围在图上呈现为一个圆盘，海口、拉萨、乌鲁木齐等地的圆盘应以对应雷达站为圆心；按 Web Mercator 解读时，这些圆盘会整体偏移到蒙古国境内。

像素格式

每块瓦片为 256x256 个大端 int16 数值，行 0 在瓦片南缘（自南向北排列）、列 0 在西缘：0x7fff（32767）表示缺测，其余数值除以 10 得到 dBZ。完全无数据的瓦片不发布，请求返回 404。请求需要携带浏览器 User-Agent，否则会被服务器拒绝。

1
2
3
4
5

def decode_tile(buf: bytes) -> np.ndarray:
    raw = np.frombuffer(buf, dtype=">i2").reshape(256, 256)
    out = raw.astype(np.float32) * 0.1
    out[raw == 32767] = np.nan
    return out[::-1]  # 翻转为北上（north-up）方向

时次发现

服务没有公开的、不依赖浏览器会话的时次列表接口，只能通过探测一个固定的参考瓦片来判断某个时次是否已发布：

1

PROBE_TILE = (3, 6, 2)  # z=3 时覆盖东经 90-135、北纬 0-45，覆盖中国大部

对该瓦片发 HEAD 请求，200 视为该时次已发布，404 视为未发布。latest_time 从当前时刻向前逐档探测，available_times 则并发探测一个时间区间内的所有档位。

缺档回填

某些时次在某些层级只发布了部分瓦片，例如 z=5 时北纬 32 度以南的行经常缺失，而 z=4 在同一区域有数据。对拼图中每个缺失的瓦片，实现里沿层级向上查找最近发布的上级瓦片（z-1、z-2……），取出上级瓦片中对应的象限区域，再用最近邻方式放大填回：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

def _fill_from_ancestors(...):
    for x, y in missing:
        for dz in range(1, zoom - MIN_ZOOM + 1):
            az, ax, ay = zoom - dz, x >> dz, y >> dz
            ancestor = get_ancestor(az, ax, ay)
            if ancestor is None:
                continue
            sub = TILE_SIZE >> dz
            row = ((1 << dz) - 1 - (y - (ay << dz))) * sub
            col = (x - (ax << dz)) * sub
            quad = ancestor[row : row + sub, col : col + sub]
            filled[(x, y)] = np.repeat(np.repeat(quad, 2**dz, axis=0), 2**dz, axis=1)
            break

row 的计算需要额外处理坐标系差异：上级瓦片解码后是北上方向，而瓦片 y 索引是自南向北计数，两者的行序相反，(1 << dz) - 1 - (y - (ay << dz)) 是这个翻转的换算。

拼接为 NetCDF

fetch_dataset 按 bbox 和 zoom 计算所需的瓦片范围，并发抓取多个时次的瓦片、拼接为二维数组，再组装为 CF-1.8 风格的 xarray.Dataset：变量 cref(time, lat, lon)，单位 dBZ；time 为 UTC，lat/lon 为瓦片像素中心经纬度。写入 NetCDF 时以 int16 加 scale_factor=0.1 打包并启用 zlib 压缩：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

encoding = {
    "cref": {
        "dtype": "int16",
        "scale_factor": 0.1,
        "_FillValue": np.int16(32767),
        "zlib": True,
        "complevel": 4,
    },
    "time": {"units": "seconds since 1970-01-01T00:00:00Z"},
}

结果

下图是 2026-08-16 04:36/04:42/04:48 UTC 三个连续时次的 zoom=4 拼图，叠加国境线与雷达覆盖圆盘（灰色区域）：

三个连续时次的拼图结果，叠加国境线与雷达覆盖圆盘

三个连续时次的拼图结果，叠加国境线与雷达覆盖圆盘

彩色像素为解码后的 dBZ 值，均落在灰色覆盖圆盘范围内，且圆盘以雷达站为圆心，与前述验证方法一致；三个时次之间的回波位置连续演变，与降水系统的实际移动相符。