使用salem处理wrfout数据，进行切片、并插值到等压面

使用salem处理wrfout数据，进行切片、并插值到等压面

最近，有学习到通过salem处理wrfout数据，非常的简单快捷，读取的变量也比较方面。也可以快速出图，下面简单对比一下xarray和salem读取wrfout文件的区别：

import xarray as xr

from salem.utils import get_demo_file

ds = xr.open_dataset(get_demo_file('wrfout_d01.nc'))

可以发现，通过xarray 读取wrfout文件的内容非常的混乱，让人摸不着头脑，无从下手

Out[4]: ds
<xarray.Dataset>
Dimensions:  (Time: 3, south_north: 150, west_east: 150, bottom_top: 27, west_east_stag: 151, south_north_stag: 151, bottom_top_stag: 28)
Coordinates:
    XLONG    (Time, south_north, west_east) float32 ...
    XLONG_U  (Time, south_north, west_east_stag) float32 ...
    XLAT_U   (Time, south_north, west_east_stag) float32 ...
    XLAT_V   (Time, south_north_stag, west_east) float32 ...
    XLONG_V  (Time, south_north_stag, west_east) float32 ...
    XLAT     (Time, south_north, west_east) float32 ...
Dimensions without coordinates: Time, south_north, west_east, bottom_top, west_east_stag, south_north_stag, bottom_top_stag
Data variables:
    Times    (Time) |S19 ...
    T2       (Time, south_north, west_east) float32 ...
    RAINC    (Time, south_north, west_east) float32 ...
    RAINNC   (Time, south_north, west_east) float32 ...
    U        (Time, bottom_top, south_north, west_east_stag) float32 ...
    V        (Time, bottom_top, south_north_stag, west_east) float32 ...
    PH       (Time, bottom_top_stag, south_north, west_east) float32 ...
    PHB      (Time, bottom_top_stag, south_north, west_east) float32 ...
Attributes:
    note:     Global attrs removed.

而如果使用salem的话，则大大简化：

import salem

ds = salem.open_wrf_dataset(get_demo_file('wrfout_d01.nc'))

Out[7]: ds
<xarray.Dataset>
Dimensions:       (time: 3, south_north: 150, west_east: 150, bottom_top: 27)
Coordinates:
    lon           (south_north, west_east) float32 70.72 70.97 ... 117.5 117.8
    lat           (south_north, west_east) float32 7.789 7.829 ... 46.52 46.46
  * time          (time) datetime64[ns] 2008-10-26T12:00:00 ... 2008-10-26T18...
  * west_east     (west_east) float64 -2.235e+06 -2.205e+06 ... 2.235e+06
  * south_north   (south_north) float64 -2.235e+06 -2.205e+06 ... 2.235e+06
Dimensions without coordinates: bottom_top
Data variables: (12/14)
    T2            (time, south_north, west_east) float32 ...
    RAINC         (time, south_north, west_east) float32 ...
    RAINNC        (time, south_north, west_east) float32 ...
    U             (time, bottom_top, south_north, west_east) float32 ...
    V             (time, bottom_top, south_north, west_east) float32 ...
    PH            (time, bottom_top, south_north, west_east) float32 ...
    ...            ...
    PRCP          (time, south_north, west_east) float32 ...
    WS            (time, bottom_top, south_north, west_east) float32 ...
    GEOPOTENTIAL  (time, bottom_top, south_north, west_east) float32 ...
    Z             (time, bottom_top, south_north, west_east) float32 ...
    PRCP_NC       (time, south_north, west_east) float32 ...
    PRCP_C        (time, south_north, west_east) float32 ...
Attributes:
    note:        Global attrs removed.
    pyproj_srs:  +proj=lcc +lat_0=29.0399971008301 +lon_0=89.8000030517578 +l...

可以发现，已经重命名了一些维度/坐标，定义了新的变量。

同时，salem也可以实现快速出图的功能:

ds.PRCP.isel(time=-1).salem.quick_map(cmap='Purples', vmax=5)

功能还是非常强大的，但是，注意的是：对于诊断量的计算并不完善，期待后续的跟进。

下面根据我的需求，使用的函数进行简单记录：

切片

首先是实现对于纬度切片的功能， 假设我只需要纬度：0-15°N范围的数据，可以这样实现：

ds = ds.salem.subset(corners=((100, 0.0), (210, 15)))

其中，corners=((100, 0.0), (210, 15))分布表示，提取矩形区域内的左下角的点以及右上角的点的经纬度坐标，如下图所示

经过我的测试，提取还是比较成功的，但是可能提取的结果并没有那么令人满意，总的来说问题不大。

插值到等压面

一般来说，wrf模式跑出的数据是： eta-levels，对于分析和作图来说是非常不方便的，因此salem也提供了 wrf_zlevel()、wrf_plevel() 两个函数，分布实现对于等高面和等压面的插值，这里仅示范插值到等压面的用法，使用如下所示：

DatasetAccessor.wrf_plevel(varname, levels=None, fill_value=nan, use_multiprocessing=True)
ds.isel(time=0).salem.wrf_plevel('V', levels=[850.])

“U”：你要插值的变量
levels：你要插值的气压层

通过上述两个方法，就能实现切片和插值的需求啦~

官网链接：salem