开源气象雷达NEXRAD数据完整获取与Python实操教程【完整版】

KTLX20260421_000653_V06

pip install requests boto3 arm-pyart matplotlib numpy cartopy

import osimport boto3import requestsfrom botocore import UNSIGNEDfrom botocore.config import Configfrom datetime import datetime, timedelta# ===================== 配置参数 =====================RADAR_ID = "KTLX"START_DT = datetime(2026, 4, 21, 12, 0, 0)END_DT   = datetime(2026, 4, 21, 14, 0, 0)SAVE_DIR = "./nexrad_level2_data"# ====================================================# 匿名访问S3（无需任何凭证）s3 = boto3.client("s3", region_name="us-east-1",                  config=Config(signature_version=UNSIGNED))BUCKET_NAME = "unidata-nexrad-level2"def get_nexrad_file_list(radar_id, target_date):    """获取指定雷达+日期下的全部Level2文件（使用正确前缀）"""    date_str = target_date.strftime("%Y/%m/%d")    prefix = f"{date_str}/{radar_id}/"   # 关键：日期在前，雷达ID在后    file_list = []    resp = s3.list_objects_v2(Bucket=BUCKET_NAME, Prefix=prefix)    if "Contents" not in resp:        return file_list    for obj in resp["Contents"]:        key = obj["Key"]        if key.endswith(("V06", "V07")):            file_list.append(key)    return file_listdef parse_file_time(file_key):    """从文件名解析观测时间（适配KTLX20260421_000653_V06格式）"""    fname = os.path.basename(file_key)          # 如 KTLX20260421_000653_V06    parts = fname.split('_')    if len(parts) < 2:        return None    # 第一部分如 KTLX20260421，取后8位日期    date_part = parts[0][-8:]                   # 20260421    time_part = parts[1]                        # 000653    try:        file_dt = datetime.strptime(date_part + time_part, "%Y%m%d%H%M%S")        return file_dt    except Exception:        return Nonedef download_s3_file(file_key, save_path):    """流式下载S3文件"""    url = f"https://{BUCKET_NAME}.s3.amazonaws.com/{file_key}"    with requests.get(url, stream=True, timeout=60) as r:        r.raise_for_status()        with open(save_path, "wb") as f:            for chunk in r.iter_content(chunk_size=1024*1024):                f.write(chunk)    print(f"✅ 下载完成: {os.path.basename(save_path)}")def main():    os.makedirs(SAVE_DIR, exist_ok=True)    day_range = (END_DT.date() - START_DT.date()).days + 1    for delta_day in range(day_range):        cur_date = START_DT.date() + timedelta(days=delta_day)        print(f"正在扫描日期: {cur_date}")        file_keys = get_nexrad_file_list(RADAR_ID, cur_date)        print(f"  找到 {len(file_keys)} 个候选文件")        for key in file_keys:            f_dt = parse_file_time(key)            if f_dt is None:                continue            if START_DT <= f_dt <= END_DT:                fname = os.path.basename(key)                save_path = os.path.join(SAVE_DIR, fname)                if os.path.exists(save_path):                    print(f"⏭️ 文件已存在: {fname}")                    continue                download_s3_file(key, save_path)if __name__ == "__main__":    main()

import osimport pyartimport matplotlib.pyplot as plt# ===================== 配置 =====================DATA_DIR = "./nexrad_level2_data"# ================================================def plot_radar_ppi(file_path):    # 读取NEXRAD Level 2数据（直接读取，pyart自动处理压缩）    radar = pyart.io.read(file_path)    print("=" * 50)    print("雷达站点：", radar.metadata["instrument_name"])    print("观测时间：", radar.time["units"])    print("可用观测字段：", list(radar.fields.keys()))    print("扫描仰角数量：", radar.nsweeps)    # 选取第一个仰角进行绘图    sweep_idx = 0    sweep_angle = radar.fixed_angle['data'][sweep_idx]    # 创建画布    plt.figure(figsize=(10, 9), dpi=120)    display = pyart.graph.RadarDisplay(radar)    # 绘制反射率PPI图    display.plot(        field="reflectivity",        sweep=sweep_idx,        vmin=-10, vmax=70,        cmap="NWSRef",        title=f"NEXRAD {radar.metadata['instrument_name']}\n"              f"Reflectivity PPI | Elevation {sweep_angle:.2f}°"    )    # 设置绘图范围为150km    display.set_limits(xlim=(-150, 150), ylim=(-150, 150))    plt.tight_layout()    plt.show()def main():    file_list = [f for f in os.listdir(DATA_DIR) if not os.path.isdir(os.path.join(DATA_DIR, f))]    # 注意：下载的文件无扩展名，直接按文件名过滤即可    if not file_list:        print("❌ 未找到雷达数据文件，请先运行下载代码")        return    first_file = os.path.join(DATA_DIR, file_list[0])    plot_radar_ppi(first_file)if __name__ == "__main__":    main()

import osimport pyartimport matplotlib.pyplot as pltimport numpy as np# ================== 配置参数 ==================DATA_DIR = "./nexrad_level2_data"      # 数据目录SWEEP_IDX = 0                          # 仰角索引（0为最低仰角）RANGE_KM = 150                         # 显示半径（公里）OUTPUT_DIR = "."                       # 图片保存目录VALID_THRESHOLD = 5.0                  # 有效数据占比阈值（%），低于此值跳过绘图# ============================================def check_field_validity(radar, field_name, sweep_idx):    """    检查指定字段在指定仰角下的有效数据占比    返回: (valid_ratio, data_slice)    """    if field_name not in radar.fields:        return 0.0, None    sweep_slice = radar.get_slice(sweep_idx)    data = radar.fields[field_name]['data'][sweep_slice]    total = data.size    if total == 0:        return 0.0, None    if hasattr(data, 'mask'):        valid = (~data.mask).sum()    else:        valid = total    ratio = 100.0 * valid / total    return ratio, sweep_slicedef plot_multi_field_ppi(file_path, sweep_idx=0, range_km=150, output_dir="."):    """    绘制多参数 PPI 对比图并保存    返回: bool, 是否成功生成图片    """    fname = os.path.basename(file_path)    try:        radar = pyart.io.read(file_path)    except Exception as e:        print(f"  ⚠️ 读取失败，跳过文件: {fname} | 错误: {str(e)[:80]}")        return False    # 站点信息和仰角    site_name = radar.metadata.get('instrument_name', 'Unknown')    if sweep_idx >= radar.nsweeps:        print(f"  ⚠️ 仰角索引 {sweep_idx} 超出范围 (共{radar.nsweeps}个)，使用仰角0")        sweep_idx = 0    sweep_angle = radar.fixed_angle['data'][sweep_idx]    # 定义想要绘制的字段（英文标签）    fields_config = {        'reflectivity': {            'label': 'Reflectivity (dBZ)',            'vmin': -20, 'vmax': 70,            'cmap': 'NWSRef',        },        'velocity': {            'label': 'Radial Velocity (m/s)',            'vmin': -30, 'vmax': 30,            'cmap': 'NWSVel',        },        'spectrum_width': {            'label': 'Spectrum Width (m/s)',            'vmin': 0, 'vmax': 10,            'cmap': 'NWS_SPW',        },        'differential_reflectivity': {            'label': 'Differential Reflectivity ZDR (dB)',            'vmin': -2, 'vmax': 6,            'cmap': 'RefDiff',        }    }    # 筛选出文件中存在且有效数据占比大于阈值的字段    available_fields = []    for field in fields_config:        ratio, _ = check_field_validity(radar, field, sweep_idx)        if ratio >= VALID_THRESHOLD:            available_fields.append(field)        elif field in radar.fields:            print(f"  ⚠️ 字段 {field} 有效数据占比 {ratio:.1f}% < {VALID_THRESHOLD}%，跳过绘制")        else:            print(f"  ⚠️ 字段 {field} 不存在，跳过")    if not available_fields:        print(f"  ❌ 无有效字段可绘制，跳过 {fname}")        return False    n = len(available_fields)    fig, axes = plt.subplots(1, n, figsize=(5 * n, 5))    if n == 1:        axes = [axes]    for ax, field_name in zip(axes, available_fields):        cfg = fields_config[field_name]        display = pyart.graph.RadarDisplay(radar)        try:            # 绘制 PPI            display.plot_ppi(                field_name, sweep=sweep_idx,                vmin=cfg['vmin'], vmax=cfg['vmax'],                cmap=cfg['cmap'], ax=ax,                colorbar_flag=True,                title_flag=False            )            # 设置显示范围            display.set_limits(xlim=(-range_km, range_km), ylim=(-range_km, range_km), ax=ax)            ax.set_title(cfg['label'], fontsize=12)            # 添加仰角信息            ax.text(0.02, 0.95, f"Elev {sweep_angle:.1f}°",                    transform=ax.transAxes, fontsize=9,                    verticalalignment='top', bbox=dict(facecolor='white', alpha=0.7))        except Exception as e:            print(f"  ⚠️ 绘制 {field_name} 失败: {e}")            ax.text(0.5, 0.5, f"{field_name}\nplot error", ha='center', transform=ax.transAxes)            ax.set_title(cfg['label'], fontsize=12)    # 总标题    fig.suptitle(f"NEXRAD {site_name} | Multi-Parameter PPI (sweep {sweep_idx})", fontsize=14, y=1.02)    plt.tight_layout()    # 保存图片    base = os.path.splitext(fname)[0]    out_filename = f"{base}_sweep{sweep_idx:02d}_multi.png"    out_path = os.path.join(output_dir, out_filename)    plt.savefig(out_path, dpi=150, bbox_inches='tight')    print(f"  ✅ 已保存: {out_path}")    plt.close(fig)    return Truedef main():    if not os.path.exists(DATA_DIR):        print(f"错误：数据目录 {DATA_DIR} 不存在，请先下载数据。")        return    files = [f for f in os.listdir(DATA_DIR)             if os.path.isfile(os.path.join(DATA_DIR, f))]    if not files:        print(f"错误：数据目录 {DATA_DIR} 为空。")        return    print(f"找到 {len(files)} 个文件，开始处理...\n")    for idx, fname in enumerate(files, 1):        print(f"[{idx}/{len(files)}] 处理: {fname}")        full_path = os.path.join(DATA_DIR, fname)        plot_multi_field_ppi(full_path, sweep_idx=SWEEP_IDX,                             range_km=RANGE_KM, output_dir=OUTPUT_DIR)        print()    print("所有文件处理完毕。")if __name__ == "__main__":    main()