我正在尝试将Lambert共形数据集重新投影到Plate Carree。我知道这可以很容易地使用卡托皮直观地完成。但是,我正在尝试创建一个新数据集,而不仅仅是显示重新投影的图像。以下是我绘制的方法,但我无法正确子集数据集(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 插值从t850_subset创建一个新的平方纬度/纬度数据集。
有没有更简单的方法来调整 xarray 数据集的大小和重新投影?
你最简单的方法就是利用xarray的能力来做类似熊猫的数据选择;这是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还有一套其他插值方法。
我们在 Iris(鸢尾花)中有多种"重新网格化"方法,如果这对您来说不是太多的上下文切换。
Xarray 在这里解释了它与 Iris 的关系,并提供了一个 to_iris() 方法。