2

我正在尝试将 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 数据集？

4

2 回答 2

3

您最简单的前进道路是利用 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 的选择和索引功能的文档中找到更多信息。您可能还会对 xarray对插值的内置支持感兴趣。如果对 SciPy 之外的插值方法感兴趣，MetPy 还有一套其他插值方法。

于 2019-01-17T22:22:24.793 回答
1

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

于 2019-01-15T10:53:04.563 回答