使用append之后数组维度消失_气象编程 | 使用 Dask 并行抽取气象站点数据

本文介绍如何使用 Dask 从多个文件中提取站点垂直廓线和剖面图数据。

简介

最近在构建面向冬奥会服务的 GRAPES MESO 1KM 产品后处理的 ecFlow 系统 (grapes_meso_1km_post)。

模式生成数据的信息如下：

• 纬度 latitude：44.5 - 34.5

• 经度 longitude：108 - 124

• 分辨率：0.01 度，1km

• 格点数 ni * nj：1601 * 1601

• 时效 step：0 - 24 小时，逐小时间隔

• 层次 pl：15 个等压面层，从 1000 hPa 到 100 hPa

其中一种产品是站点预报产品，包括：

• 站点的垂直廓线数据

• 经过站点的沿经线和纬线的垂直剖面图数据

该数据由原始分辨率 GRIB 2 数据文件抽取生成，所有时效和所有变量都包含在同一个 NetCDF 文件中。一共有 9 种变量，每种变量 3 类数据，每类数据 25 个时效，每个时效包含 15 层。

当前任务脚本使用 NCL 抽取数据并输出 NetCDF 文件。经集成测试，每个站点的任务运行耗时均在 70 分钟左右。对于计划逐小时运行的模式系统来说，这 可能 不是一个理想的运行时间。

因为 ecFlow 系统需要保证前一天最后一个任务的完成时间要早于下一天第一个任务的启动时间，才能持续滚动循环。

当然，任何系统都有实现的方式。只要将系统拆分为多个 ecFlow 系统，就可以规避时间重叠的问题。只不过这通常不是最佳实践。

本文首先介绍 NCL 脚本的实现方法，再介绍笔者使用 Dask 库进行并行抽取的尝试。

使用 NCL 提取

GRAPES 模式输出的 GRIB 2 文件分时效保存。NCL 脚本首先使用 addfiles 命令批量打开 GRIB 2 文件。

inputf = addfiles(input_filename,"r")

为所有要素场申请空间。下面代码是位势高度的示例代码，其中 hgt_0 和 hgt_9 保存剖面图数据，hgt 保存廓线图数据。所有数据前两个维度都是时效和层次。

hgt_0=new((/dimsizes(time),dimsizes(level),interp_0/),"float")
hgt_9=new((/dimsizes(time),dimsizes(level),interp_9/),"float")
hgt=new((/dimsizes(time),dimsizes(level)/),"float")

从每个时效文件中依次读取要素场，填充到数据数组中。

do i=0,dimsizes(time)-1
    hgt_0(i,:,:)=inputf[i]->HGT_P0_L100_GLL0(:,180:860:-1,797)
    hgt_9(i,:,:)=inputf[i]->HGT_P0_L100_GLL0(:,405,520:1220)
    hgt(i,:)=inputf[i]->HGT_P0_L100_GLL0(:,405