从GrADS文件中加载要素场

programming/python, meteorology/data
目录

本文介绍如何使用 nwpc-oper/nwpc-data 库从 GrADS 格点二进制数据文件中加载单个要素场。

API

nwpc_data.format.grads.load_field_from_file() 函数用于从 GrADS 格点数据文件中加载要素场。

def load_field_from_file(
        file_path: Union[str, Path],
        parameter: str = None,
        level_type: str = None,
        level: Union[int, float, List] = None,
        level_dim: Optional[str] = None,
        latitude_direction: str = "degree_north",
        forecast_time: Union[str, pd.Timedelta] = None,
        **kwargs
) -> Optional[xr.DataArray]:
    pass

参数

parameter：要素名

level_type 与 level：层次类型和层次值

  • level_type 为 pl 或 ml，表示气压层或模式层，level 为层次数值
  • level_type 为 single，表示单层变量
  • level_type 为 index，表示按层次序号索引，level 为层次序号，从 0 开始

latitude_direction：纬度方向，默认从北到南 (degree_north)，与 NWPC 的 GRIB 2 数据保持一致

forecast_time：预报时效，用于单个数据描述文件对应多个单时效数据文件的情况，例如 GRAPES TYM

返回值

xarray.DataArray 对象二维格点数组，与 nwpc-data 库 GRIB 2 系列 API 返回的 DataArray 对象有相同的坐标：

  • latitude：纬度
  • longitude：经度
  • level，pl 或 ml：层次
  • valid_time：预报时间
  • start_time：起报时间
  • forecast_time：预报时效

使用

使用 find_local_file() 获取 GRAPES GFS 预报系统等压面 GrADS 数据文件路径：

from nwpc_data.data_finder.local import find_local_file

postvar_file_path = find_local_file(
    "grapes_gfs_gmf/bin/postvar_ctl",
    start_time="2021092600",
    forecast_time="36h",
    storage_base="M:"
)
print(postvar_file_path)

M:\GRAPES_GFS_GMF\Fcst-long\2021092600\post.ctl_2021092600_036

等压面要素

单个要素场：850 hPa 温度场 (level=850)

from nwpc_data.format.grads import load_field_from_file

field = load_field_from_file(
    postvar_file_path,
    parameter="t",
    level_type="pl",
    level=850
)
print(field)

<xarray.DataArray 't' (latitude: 720, longitude: 1440)>
array([[262.11594, 262.09723, 262.10767, ..., 262.10422, 262.10376,
        262.10025],
       [262.0563 , 262.0463 , 262.04837, ..., 262.04083, 262.03967,
        262.05145],
       [262.0379 , 262.04916, 262.04547, ..., 262.0374 , 262.03357,
        262.03827],
       ...,
       [224.35977, 224.34584, 224.35239, ..., 224.33386, 224.33951,
        224.37468],
       [224.59709, 224.52245, 224.46587, ..., 224.68697, 224.47331,
        224.51013],
       [225.364  , 225.61134, 225.5598 , ..., 225.22418, 225.36223,
        225.35268]], dtype=float32)
Coordinates:
  * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    pl             int32 850
    valid_time     datetime64[ns] 2021-09-27T12:00:00
    start_time     datetime64[ns] 2021-09-26
    forecast_time  timedelta64[ns] 1 days 12:00:00
Attributes:
    description:  temperature

多个要素场：500 hPa 和 850 hPa 温度场 (level=[850, 500])

field = load_field_from_file(
    postvar_file_path,
    parameter="t",
    level_type="pl",
    level=[850, 500],
)
print(field)

<xarray.DataArray 't' (pl: 2, latitude: 720, longitude: 1440)>
array([[[262.11594, 262.09723, 262.10767, ..., 262.10422, 262.10376,
         262.10025],
        [262.0563 , 262.0463 , 262.04837, ..., 262.04083, 262.03967,
         262.05145],
        [262.0379 , 262.04916, 262.04547, ..., 262.0374 , 262.03357,
         262.03827],
        ...,
        [224.35977, 224.34584, 224.35239, ..., 224.33386, 224.33951,
         224.37468],
        [224.59709, 224.52245, 224.46587, ..., 224.68697, 224.47331,
         224.51013],
        [225.364  , 225.61134, 225.5598 , ..., 225.22418, 225.36223,
         225.35268]],
       [[240.21898, 240.21875, 240.20212, ..., 240.21648, 240.19972,
         240.2088 ],
        [240.27127, 240.2787 , 240.27321, ..., 240.27722, 240.26988,
         240.26624],
        [240.3347 , 240.34404, 240.34558, ..., 240.33405, 240.3301 ,
         240.32549],
        ...,
        [227.47325, 227.47318, 227.47247, ..., 227.4741 , 227.47241,
         227.47177],
        [227.49261, 227.49744, 227.4924 , ..., 227.49075, 227.49005,
         227.487  ],
        [227.5075 , 227.52118, 227.50984, ..., 227.51366, 227.49881,
         227.49263]]], dtype=float32)
Coordinates:
  * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
  * pl             (pl) int64 850 500
    valid_time     datetime64[ns] 2021-09-27T12:00:00
    start_time     datetime64[ns] 2021-09-26
    forecast_time  timedelta64[ns] 1 days 12:00:00
Attributes:
    description:  temperature

单层要素

海平面气压

field = load_field_from_file(
    postvar_file_path,
    parameter="psl",
    level_type="single",
)
print(field)

<xarray.DataArray 'psl' (latitude: 720, longitude: 1440)>
array([[1028.3733 , 1028.3723 , 1028.3717 , ..., 1028.3759 , 1028.3746 ,
        1028.3738 ],
       [1028.3737 , 1028.3713 , 1028.3694 , ..., 1028.38   , 1028.3777 ,
        1028.3757 ],
       [1028.2827 , 1028.2795 , 1028.2749 , ..., 1028.2948 , 1028.2904 ,
        1028.2866 ],
       ...,
       [ 984.6077 ,  984.6296 ,  984.6277 , ...,  984.62537,  984.6316 ,
         984.6019 ],
       [ 984.00903,  984.0521 ,  984.0917 , ...,  983.9589 ,  984.0802 ,
         984.05743],
       [ 983.55194,  983.4115 ,  983.44916, ...,  983.6367 ,  983.5626 ,
         983.54865]], dtype=float32)
Coordinates:
  * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    level          float64 0.0
    valid_time     datetime64[ns] 2021-09-27T12:00:00
    start_time     datetime64[ns] 2021-09-26
    forecast_time  timedelta64[ns] 1 days 12:00:00
Attributes:
    description:  sea level pressure

多层要素

土壤温度

注：结果中 level 值有问题，应该使用层次序号代替，有待后续修正

field = load_field_from_file(
    postvar_file_path,
    parameter="tsoil",
    level_type="index",
    level=0
)
print(field)

<xarray.DataArray 'tsoil' (latitude: 720, longitude: 1440)>
array([[255.38072, 255.39757, 255.3218 , ..., 255.43384, 255.37071,
        255.34639],
       [255.02567, 255.00705, 255.00714, ..., 255.06354, 254.96187,
        254.96585],
       [255.12096, 255.08827, 255.12408, ..., 255.09755, 255.11043,
        255.1201 ],
       ...,
       [218.36595, 217.44775, 217.25752, ..., 218.07549, 217.9659 ,
        217.98798],
       [218.30225, 217.31868, 217.17842, ..., 217.97504, 217.88512,
        217.87366],
       [218.21353, 217.29764, 217.18007, ..., 217.9085 , 217.83109,
        217.77658]], dtype=float32)
Coordinates:
  * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    level          float64 1e+03
    valid_time     datetime64[ns] 2021-09-27T12:00:00
    start_time     datetime64[ns] 2021-09-26
    forecast_time  timedelta64[ns] 1 days 12:00:00
Attributes:
    description:  soil temperature

模式层

获取 GRAPES GFS 模式层 GrADS 输出的数据描述文件路径：

modelvar_file_path = find_local_file(
    "grapes_gfs_gmf/bin/modelvar_ctl",
    start_time="2021092600",
    forecast_time="120h",
    storage_base="M:"
)
print(modelvar_file_path)

M:\GRAPES_GFS_GMF\Fcst-long\2021092600\model.ctl_2021092600_120

加载第 10 模式层的 u 分量：

field = load_field_from_file(
    modelvar_file_path,
    parameter="u",
    level_type="ml",
    level=10
)
print(field)

<xarray.DataArray 'u' (latitude: 721, longitude: 1440)>
array([[6.4011283, 6.374322 , 6.3444586, ..., 6.4519706, 6.422714 ,
        6.40848  ],
       [5.8047924, 5.7803693, 5.7510133, ..., 5.919187 , 5.8902326,
        5.8443203],
       [7.020947 , 6.9858413, 6.9516716, ..., 7.137582 , 7.1022863,
        7.060398 ],
       ...,
       [4.8194976, 4.823494 , 4.830612 , ..., 4.811972 , 4.815815 ,
        4.8179264],
       [5.3669686, 5.369134 , 5.370997 , ..., 5.350346 , 5.3530393,
        5.3615303],
       [6.2383065, 6.233473 , 6.2309933, ..., 6.2335534, 6.227803 ,
        6.231739 ]], dtype=float32)
Coordinates:
  * latitude       (latitude) float64 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    ml             int32 10
    valid_time     datetime64[ns] 2021-10-01
    start_time     datetime64[ns] 2021-09-26
    forecast_time  timedelta64[ns] 20:00:00
Attributes:
    description:  U wind component

单一描述文件对应多个数据文件

GRAPES TYM 等压面 GrADS 数据每个时次只有一个描述文件，对应多个单时效二进制数据文件。

获取 POSTVAR 文件路径：

postvar_file_path = find_local_file(
    "grapes_tym/bin/postvar_ctl",
    start_time="2021092600",
    storage_base="M:"
)

M:\GRAPES_TYM\Fcst-main\2021092600\post.ctl_2021092600

CTL 文件名只包含起报日期 (2021.09.26) 和起报时次 (00)。

加载 850 hPa 温度场：

field = load_field_from_file(
    postvar_file_path,
    parameter="t",
    level_type="pl",
    forecast_time="24h",
    level=850
)

<xarray.DataArray 't' (latitude: 835, longitude: 1557)>
array([[272.97653, 272.9475 , 272.91296, ..., 268.50092, 268.47885,
        268.45917],
       [272.9719 , 272.97913, 272.9622 , ..., 268.40283, 268.38678,
        268.36868],
       [273.06686, 273.05692, 273.03488, ..., 268.2895 , 268.27045,
        268.2581 ],
       ...,
       [291.57855, 291.59143, 291.58658, ..., 289.32324, 289.31635,
        289.30954],
       [291.59384, 291.5858 , 291.56274, ..., 289.29056, 289.29034,
        289.29037],
       [291.60318, 291.5816 , 291.55527, ..., 289.2471 , 289.28748,
        289.29752]], dtype=float32)
Coordinates:
  * latitude    (latitude) float64 60.06 59.97 59.88 ... -14.82 -14.91 -15.0
  * longitude   (longitude) float64 40.0 40.09 40.18 40.27 ... 179.9 179.9 180.0
    pl          int64 850
    valid_time  datetime64[ns] 2021-09-27
Attributes:
    description:  temperature

实现

nwpc-data 库内部包含一个简单的 GrADS 数据文件解析器，包括以下几个类：

  • GradsCtlParser：用于解析 GrADS 数据描述文件
  • GradsCtl：表示 GrADS 数据描述文件，可由 GradsCtlParser 生成
  • GradsDataHandler：用于解码 GrADS 格点二进制数据文件，需要 GradsCtl 对象
  • GradsRecordHandler：表示数据文件中的一条记录，也就是单个要素场，可由 GradsDataHandler 生成

下图是从 GrADS 格点二进制数据文件中加载单个要素场的示意图。

从 GrADS 文件中加载要素场示意图

输入

  • File Path：数据描述文件路径
  • Field Filters：要素场筛选条件，即上述 API 中 parameter，level_type，level 等参数

输出

Data Array：xarray.DataArray 二维数组对象

加载步骤

  1. 使用 GradsCtlParser 解析数据描述文件 (File Path)，生成数据描述对象 GradsCtl；
  2. 根据 GradsCtl 创建数据解码器 GradsDataHandler，按要素场筛选条件 (Field Filters) 使用 find_record() 函数查找需要的要素场，得到数据记录对象 GradsRecordHandler；
  3. 调用 create_data_array_from_record() 函数，将数据记录 GradsRecordHandler 转为 xarray.DataArray 对象并返回。

参考

nwpc-oper/nwpc-data

关于 GrADS 格点数据格式的更多信息请浏览以下两篇文章：

《GrADS格点数据格式》

《GrADS格点数据文件解析》