modelwhale 格点转格点

https://www.heywhale.com/mw/project/5f2c00d5af3980002cb2ad99/content

import glob
import math
import multiprocessing
import os

import time
import warnings
from datetime import datetime, timedelta
from math import acos, asin, cos, pi, radians, sin, sqrt

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib import colors
from matplotlib import pyplot as plt
from matplotlib.cm import get_cmap
from scipy import interpolate
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import xarray as xr
import xesmf as xe
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from netCDF4 import Dataset
from wrf import ALL_TIMES, cartopy_xlim, cartopy_ylim, get_cartopy, getvar, latlon_coords, smooth2d, to_np
import warnings
warnings.filterwarnings(‘ignore’)
%matplotlib inline

格点数据到格点数据
案例讲解：将WRF模式网格调整为等经纬度网格
中尺度天气预报模式WRF，是气象研究，特别是中尺度天气研究中最长使用的数值模式。它与前文介绍的那些预报数据和再分析数据有两个很大的区别:

其预报或模拟范围往往不是全球，而是一块区域。
其网格投影方式不是常见经纬度网格。对我国大部分地区进行预报或模拟时，使用的通常是兰博特投影。
因此，当我们希望WRF网格能与fnl等全球模式网格向匹配/对齐的时候，通常我们需要以WRF的网格投影或全球模式的网格投影，两者之一为基础，另一方通过一种名为regridding(重栅格化)的操作，转化为同一种投影。

在python中，我们常使用xESMF库进行此操作。

假设这样一个业务场景:我们已经进行了一次WRF模拟或预测，现在我们需要对其进行检验。我们可以选择读取再分析资料（如fnl),将其regridding到wrf的空间网格投影上。

首先，读取wrf模式结果。wrf模式的运行结果一般都是以wrfout为开头的nc文件

wrffile = glob.glob("/home/kesci/input/wrf5525/wrfout_d01_2019-10-01_12:00:00")[0]
ds_wrf = xr.open_dataset(wrffile)

T_wrf = ds_wrf.T

从上个cell可以看出，温度T变量存在时间、高度层（实际上是eta坐标层）、纬度(实际上更准确的说法是从南到北的网格数）、经度(类似地，从西到东的网格数)。

可以看出，wrf的变量维度和fnl、era5等再分析资料是很不相同的：首先，它的经纬度坐标是二维的。这是因为中纬度地区使用WRF进行天气预报与天气模拟时，常常使用的是兰勃特投影（这在下面可视化环节可以看出），而非经纬度网格；其次，它的纬度是单调递增的——这正好与再分析资料相反。

类似地，我们可以读取wrf的经纬度坐标变量XLAT和XLONG，并进行查看

XLAT = ds_wrf.XLAT.values[0, :, :]
XLONG = ds_wrf.XLONG.values[0, :, :]

使用xEMSF可以将两种不同网格的数据再栅格化到同一网格。

首先需要读取两种不同数据各自的经纬度网格。wrf之前已经读取过。我们来读取fnl。

fnlfile = “/home/kesci/input/fnl9943/fnl/2019/01/01/fnl_20190101_00_00.grib2”
ds_fnl = xr.open_dataset(fnlfile,engine=‘pynio’)
T_fnl = ds_fnl.TMP_P0_L100_GLL0[0,…]

接着需要做的是，将两个数据的坐标名称设置相同。可以直接通过新建xarray.Dataset的方法加以实现：

输出目标，wrf的网格

ds_out = xr.Dataset(
{
‘lon’: ([‘south_north’, ‘west_east’], XLONG),
‘lat’: ([‘south_north’, ‘west_east’], XLAT),
}
)
ds_out

输入数据，fnl网格

ds_in = xr.Dataset(
{
‘lat’: ([‘west_east’], ds_fnl.lat_0),
‘lon’: ([‘south_north’], ds_fnl.lon_0),
}
)
ds_in

然后，将两组网格作为变量输入xe.Regridder函数，使用双线性插值，生成一个regridder

regridder = xe.Regridder(ds_in, ds_out, ‘bilinear’, reuse_weights=False) #reuse_weights=False表示不保存权重系数到一个文件

然后将fnl模式的变量作为参数，输入此regridder，生成的就是一个wrf网格的新变量T_final

T_final = regridder(T_fnl)