重新投影 Xarray 数据集答案

【问题标题】：Reprojecting Xarray Dataset重新投影 Xarray 数据集
【发布时间】：2019-06-08 21:19:47
【问题描述】：

我正在尝试将 Lambert Conformal 数据集重新投影到 Plate Carree。我知道这可以很容易地使用 cartopy 在视觉上完成。但是，我正在尝试创建一个新数据集，而不仅仅是显示重新投影的图像。以下是我制定的方法，但我无法正确地对数据集进行子集化（Python 3.5、MacOSx）。

from siphon.catalog import TDSCatalog
import xarray as xr
from xarray.backends import NetCDF4DataStore
import numpy as np
import cartopy.crs as ccrs
from scipy.interpolate import griddata
import numpy.ma as ma
from pyproj import Proj, transform
import metpy

# Declare bounding box
min_lon = -78
min_lat = 36
max_lat = 40
max_lon = -72
boundinglat = [min_lat, max_lat]
boundinglon = [min_lon, max_lon]

# Load the dataset
cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/HRRR/CONUS_2p5km/latest.xml')
dataset_name = sorted(cat.datasets.keys())[-1]
dataset = cat.datasets[dataset_name]
ds = dataset.remote_access(service='OPENDAP')
ds = NetCDF4DataStore(ds)
ds = xr.open_dataset(ds)

# parse the temperature at 850 and @ 0z reftime
tempiso = ds.metpy.parse_cf('Temperature_isobaric')
t850 = tempiso[0][2]

# transform bounding lat/lons to src_proj
src_proj = tempiso.metpy.cartopy_crs #aka lambert conformal conical
extents = src_proj.transform_points(ccrs.PlateCarree(), np.array(boundinglon), np.array(boundinglat))

# subset the data using the indexes of the closest values to the src_proj extents
t850_subset = t850[(np.abs(tempiso.y.values - extents[1][0])).argmin():(np.abs(tempiso.y.values - extents[1][1])).argmin()][(np.abs(tempiso.x.values - extents[0][1])).argmin():(np.abs(tempiso.x.values - extents[0][0])).argmin()]

# t850_subset should be a small, reshaped dataset, but it's shape is 0x2145
# now use nplinspace, npmeshgrid & scipy interpolate to reproject

我的变换点 > 查找最近值子集不起作用。它声称最近的点在数据集的范围之外。如前所述，我计划使用 nplinspace、npmeshgrid 和 scipy interpolate 从 t850_subset 创建一个新的方形 lat/lon 数据集。

有没有更简单的方法来调整和重新投影 xarray 数据集？

【问题讨论】：

获取的数据集的形状（numrows, numcols）是什么？
@swatchai 形状为 1377 x 2145
如果可以将t850转换为numpy的数组，那么可以对数组进行裁剪或切片得到子集部分。

标签： python-3.x python-xarray cartopy pyproj metpy

【解决方案1】：

您最简单的方法是利用 xarray 进行类似 pandas 的数据选择的能力；这是 IMO 最好的 xarray 部分。将最后两行替换为：

# By transposing the result of transform_points, we can unpack the
# x and y coordinates into individual arrays.
x_lim, y_lim, _ = src_proj.transform_points(ccrs.PlateCarree(),
    np.array(boundinglon), np.array(boundinglat)).T
t850_subset = t850.sel(x=slice(*x_lim), y=slice(*y_lim))

您可以在有关 xarray 的 selection and indexing functionality 的文档中找到更多信息。您可能还会对 xarray 的built-in support for interpolation 感兴趣。如果对 SciPy 之外的插值方法感兴趣，MetPy 也有一套 other interpolation methods。

【讨论】：

【解决方案2】：

Iris 中有各种“重新网格化”方法，如果这对您来说不是太多的上下文切换。
Xarray 解释了它与 Iris here 的关系，并提供了a to_iris() method。

【讨论】：