重新投影 Xarray 数据集
Reprojecting Xarray Dataset
我正在尝试将 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 插值从 t850_subset 创建一个新的正方形 lat/lon 数据集。
是否有更简单的方法来调整和重新投影 xarray 数据集?
我们在 Iris 中提供了多种 "regridding" 方法,如果这对您来说不是太多的上下文切换的话。
Xarray 解释它与 Iris 的关系 here, and provides a to_iris() method.
您最简单的前进道路是利用 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. You would probably also be interested in xarray's built-in support for interpolation. And if interpolation methods beyond SciPy's are of interest, MetPy also has a suite of other interpolation methods.
文档中找到更多信息
我正在尝试将 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 插值从 t850_subset 创建一个新的正方形 lat/lon 数据集。
是否有更简单的方法来调整和重新投影 xarray 数据集?
我们在 Iris 中提供了多种 "regridding" 方法,如果这对您来说不是太多的上下文切换的话。
Xarray 解释它与 Iris 的关系 here, and provides a to_iris() method.
您最简单的前进道路是利用 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. You would probably also be interested in xarray's built-in support for interpolation. And if interpolation methods beyond SciPy's are of interest, MetPy also has a suite of other interpolation methods.
文档中找到更多信息