xwrf | 批量处理多个wrf文件并插值指定高度层

用户11172986

发布于 2026-04-24 18:40:13

610

xwrf | 批量处理多个wrf文件并插值指定高度层

项目概述

之前有发出节目预告，不知各位有没看到公众号的文章。

WRF（Weather Research and Forecasting）模式是广泛应用于气象研究和业务预报的中尺度数值天气预报系统，其输出结果通常包含多维气象要素场。在实际科研和业务应用中，我们经常需要处理多个WRF输出文件（通常以NetCDF格式存储），并对这些数据进行后处理和分析。其中，将三维气象要素插值到指定高度层是一项常见且重要的需求，这有助于在不同高度层次上分析大气状态，便于垂直剖面分析或与其他观测数据进行对比。

本教程将介绍如何使用xWRF工具（基于xarray的WRF数据处理扩展）批量处理多个WRF输出文件，并将三维气象场插值到用户指定的高度层。这种方法相比传统逐文件处理方式更加高效，能够充分利用Python生态中的并行处理能力，同时保持代码的简洁性和可读性。

环境设置

安装依赖

!pip install xwrf -i https://pypi.mirrors.ustc.edu.cn/simple/

!pip install pint_xarray  -i https://pypi.mirrors.ustc.edu.cn/simple/

!pip install xgcm -i https://pypi.mirrors.ustc.edu.cn/simple/

导入必要的库

import xarray as xr
import xwrf
# for unit conversion (`pint` DataArray accessor)
import pint_xarray
# for `dask`-accelerated computation
from distributed import LocalCluster, Client
import numpy as np
import glob

数据处理

files = sorted(glob.glob('/home/mw/input/typhoon9537/*'))

wrf_ds = xr.open_mfdataset(
    files,
    parallel=True,
    concat_dim="Time",
    combine="nested",
    chunks={'Time': 1},
)
print(wrf_ds)

import time
start = time.time()
xwrf_ds = wrf_ds.xwrf.postprocess()
end = time.time()
print(f"{end - start} s elapsed")
print(xwrf_ds)

0.21268272399902344 s elapsed

基于dask对于批量读取的wrf数据进行坐标转换，速度很快

交错网格处理（风变量）

wrf数据是交错网格，这个我们在

本项目来源于和鲸社区，使用转载需要标注来源作者: 酷炫用户名来源: https://www.heywhale.com/mw/project/67093d064f507a297c3142f3 WRF | 为什么wrfout中经向风和纬向风的格点数不一样

中讲过，不再赘述

在wrf-python中，我们可以直接读取ua和va解决

在xwrf中则需要另外处理

print(xwrf_ds['U'])

<xarray.DataArray 'U' (Time: 13, z: 34, y: 437, x_stag: 448)> Size: 346MB
dask.array<concatenate, shape=(13, 34, 437, 448), dtype=float32, chunksize=(1, 34, 437, 448), chunktype=numpy.ndarray>
Coordinates:
    XTIME    (Time) datetime64[ns] 104B dask.array<chunksize=(1,), meta=np.ndarray>
    XLAT_U   (y, x_stag) float32 783kB dask.array<chunksize=(437, 448), meta=np.ndarray>
    XLONG_U  (y, x_stag) float32 783kB dask.array<chunksize=(437, 448), meta=np.ndarray>
  * z        (z) float32 136B 0.9965 0.988 0.9765 0.962 ... 0.033 0.0195 0.0065
  * Time     (Time) datetime64[ns] 104B 2019-08-08T18:00:00 ... 2019-08-09T06...
  * y        (y) float64 3kB -6.54e+05 -6.51e+05 -6.48e+05 ... 6.51e+05 6.54e+05
  * x_stag   (x_stag) float64 4kB -6.705e+05 -6.675e+05 ... 6.675e+05 6.705e+05
Attributes:
    FieldType:      104
    MemoryOrder:    XYZ
    description:    x-wind component
    units:          m s-1
    stagger:        X
    standard_name:  x_wind
    grid_mapping:   wrf_projection

xwrf_ds['UA'] = xwrf_ds['U'].xwrf.destagger().variable
print(xwrf_ds['UA'])

<xarray.DataArray 'UA' (Time: 13, z: 34, y: 437, x: 447)> Size: 345MB
dask.array<mul, shape=(13, 34, 437, 447), dtype=float32, chunksize=(1, 34, 437, 447), chunktype=numpy.ndarray>
Coordinates:
    XLAT     (y, x) float32 781kB dask.array<chunksize=(437, 447), meta=np.ndarray>
    XLONG    (y, x) float32 781kB dask.array<chunksize=(437, 447), meta=np.ndarray>
    XTIME    (Time) datetime64[ns] 104B dask.array<chunksize=(1,), meta=np.ndarray>
    CLAT     (y, x) float32 781kB dask.array<chunksize=(437, 447), meta=np.ndarray>
  * z        (z) float32 136B 0.9965 0.988 0.9765 0.962 ... 0.033 0.0195 0.0065
  * Time     (Time) datetime64[ns] 104B 2019-08-08T18:00:00 ... 2019-08-09T06...
  * y        (y) float64 3kB -6.54e+05 -6.51e+05 -6.48e+05 ... 6.51e+05 6.54e+05
  * x        (x) float64 4kB -6.69e+05 -6.66e+05 -6.63e+05 ... 6.66e+05 6.69e+05
Attributes:
    FieldType:      104
    MemoryOrder:    XYZ
    description:    x-wind component
    units:          m s-1
    standard_name:  x_wind
    grid_mapping:   wrf_projection

可见格点数变化了

垂直层插值

对应wrf-python的interplevel函数的垂直层插值功能

xwrf则是选择了xgcm库作为外援进行垂直层插值

处理交错网格

from metpy.calc import wind_speed
import xgcm
import numpy as np
for field in ['U', 'V']:
    xwrf_ds[field] = xwrf_ds[field].xwrf.destagger().variable
xwrf_ds = xwrf_ds.metpy.quantify()

单层插值

import xgcm
import numpy as np
import pint_xarray

target_levels = np.array([250.]) # in hPa
#air_pressure = xwrf_ds.air_pressure.pint.to('hPa').metpy.dequantify()
air_pressure = xwrf_ds.air_pressure.metpy.convert_units('hPa').metpy.dequantify()
xwrf_ds['wind_speed'] =wind_speed(xwrf_ds.U, xwrf_ds.V).metpy.dequantify()

grid = xgcm.Grid(xwrf_ds, periodic=False)

_wind_speed = grid.transform(xwrf_ds.wind_speed.metpy.dequantify(), 'Z', target_levels, target_data=air_pressure, method='log')

_wind_speed = _wind_speed.compute()

/opt/conda/lib/python3.11/site-packages/xgcm/grid.py:989:

print(_wind_speed)

<xarray.DataArray 'wind_speed' (Time: 13, y: 437, x: 447, air_pressure: 1)> Size: 20MB
array([[[[ 2.15332438],
         [ 2.05429516],
         [ 1.96162581],
         ...,
         [ 8.61106545],
         [ 8.51674219],
         [ 8.42776731]],

        [[ 2.12816252],
         [ 2.02488963],
         [ 1.92792811],
         ...,
         [ 8.91645327],
         [ 8.82431303],
         [ 8.73533366]],

        [[ 2.11298364],
         [ 2.00608941],
         [ 1.90560238],
         ...,
...
         ...,
         [15.60896739],
         [16.0766267 ],
         [15.59409096]],

        [[ 4.00628541],
         [ 4.02528941],
         [ 3.97536056],
         ...,
         [15.96039649],
         [16.50404221],
         [15.90062312]],

        [[ 4.04618942],
         [ 3.99713923],
         [ 3.95577495],
         ...,
         [15.36883414],
         [15.29024682],
         [15.2005086 ]]]])
Coordinates:
    XTIME         (Time) datetime64[ns] 104B 2019-08-08T18:00:00 ... 2019-08-...
  * Time          (Time) datetime64[ns] 104B 2019-08-08T18:00:00 ... 2019-08-...
    XLAT          (y, x) float32 781kB 20.89 20.89 20.89 ... 32.78 32.78 32.78
    XLONG         (y, x) float32 781kB 116.5 116.6 116.6 ... 130.2 130.2 130.3
    CLAT          (y, x) float32 781kB 20.89 20.89 20.89 ... 32.78 32.78 32.78
  * y             (y) float64 3kB -6.54e+05 -6.51e+05 ... 6.51e+05 6.54e+05
  * x             (x) float64 4kB -6.69e+05 -6.66e+05 ... 6.66e+05 6.69e+05
  * air_pressure  (air_pressure) float64 8B 250.0

多层插值

from metpy.calc import wind_speed
import xgcm
import numpy as np
import xarray as xr

files = sorted(glob.glob('/home/mw/input/typhoon9537/*'))

xwrf_ds = xr.open_mfdataset(
    files,
    parallel=True,
    concat_dim="Time",
    combine="nested",
    chunks={'Time': 1},
).xwrf.postprocess()

def calculate_wind_speed(dataset, target_levels):
    """计算并插值风速到指定气压层"""
    # 解交错网格并添加单位
    for field in ['U', 'V']:
        dataset[field] = dataset[field].xwrf.destagger().variable
    
    # 计算风速并确保是去量纲化的DataArray
    dataset['wind_speed'] = wind_speed(dataset.U, dataset.V).metpy.dequantify()
    print(dataset.air_pressure.pint)
    # 处理气压数据，确保是去量纲化的DataArray
#   air_pressure = dataset.air_pressure.pint.to('hPa').metpy.dequantify()
    air_pressure = dataset.air_pressure.metpy.convert_units('hPa').metpy.dequantify()
    # 创建xgcm网格
    grid = xgcm.Grid(dataset, periodic=False)
    
    # 执行插值
    result = grid.transform(
        dataset.wind_speed,
        'Z',
        target_levels,
        target_data=air_pressure,
        method='log'
    ).compute()  # 确保计算结果
    
    return result

# 使用示例
target_levels = np.array([850,725,500,250.,200])  # hPa
wind_speed_multi = calculate_wind_speed(xwrf_ds, target_levels)

<pint_xarray.accessors.PintDataArrayAccessor object at 0x7f2a16aa3510

print(wind_speed_multi)

<xarray.DataArray 'wind_speed' (Time: 13, y: 437, x: 447, air_pressure: 5)> Size: 102MB
array([[[[ 9.25609613,  8.34912411,  9.4607662 ,  2.15332438,
           3.45706554],
         [ 9.3940668 ,  8.52057594,  9.56193775,  2.05429516,
           3.56831646],
         [ 9.54101839,  8.70075069,  9.65882567,  1.96162581,
           3.69223884],
         ...,
         [22.81890334, 17.9803955 , 23.20878212,  8.61106545,
           7.77787846],
         [22.70740632, 17.86505295, 23.15872668,  8.51674219,
           7.60701875],
         [22.59427188, 17.75521747, 23.09905585,  8.42776731,
           7.4744698 ]],

        [[ 9.22080649,  8.39445186,  9.43068502,  2.12816252,
           3.34152052],
         [ 9.36072463,  8.56623378,  9.53015887,  2.02488963,
           3.45191456],
         [ 9.51037999,  8.74737881,  9.62507364,  1.92792811,
           3.57544112],
...
         [ 4.92995826,  5.77725726,  3.21912271, 15.96039649,
          20.74913331],
         [ 4.71819977,  5.79624439,  4.17212933, 16.50404221,
          21.64187518],
         [ 4.56275508,  6.07743758,  3.95430655, 15.90062312,
          21.35588073]],

        [[ 5.04970797,  3.97664761,  3.14422791,  4.04618942,
           8.55706546],
         [ 5.0531973 ,  4.00633996,  3.21976034,  3.99713923,
           8.47484533],
         [ 5.0878566 ,  4.01714974,  3.33147043,  3.95577495,
           8.39265859],
         ...,
         [ 5.47635917,  6.58052366,  3.33776199, 15.36883414,
          20.65481604],
         [ 5.50332557,  6.69873229,  3.81604314, 15.29024682,
          20.69798193],
         [ 5.1656812 ,  6.6051779 ,  3.81690317, 15.2005086 ,
          20.73190761]]]])
Coordinates:
    XLAT          (y, x) float32 781kB 20.89 20.89 20.89 ... 32.78 32.78 32.78
    XLONG         (y, x) float32 781kB 116.5 116.6 116.6 ... 130.2 130.2 130.3
    XTIME         (Time) datetime64[ns] 104B 2019-08-08T18:00:00 ... 2019-08-...
    CLAT          (y, x) float32 781kB 20.89 20.89 20.89 ... 32.78 32.78 32.78
  * Time          (Time) datetime64[ns] 104B 2019-08-08T18:00:00 ... 2019-08-...
  * y             (y) float64 3kB -6.54e+05 -6.51e+05 ... 6.51e+05 6.54e+05
  * x             (x) float64 4kB -6.69e+05 -6.66e+05 ... 6.66e+05 6.69e+05
  * air_pressure  (air_pressure) float64 40B 850.0 725.0 500.0 250.0 200.0