无法将y和x坐标分配给xarray.DataArray

Question

在完成了MetPy横截面示例之后，我试图将该示例推广到NCEP NAM-12km GRIB2文件，但没有成功。通过将我的文件的DataArray与示例文件(netCDF文件)进行比较，我发现xarray.open_dataset()没有生成x坐标和y坐标。它也没有在输入的GRIB文件中找到足够的信息来生成metpy_crs坐标。我尝试将GRIB文件转换为netCDF，但这并没有解决问题。

即使将输入文件视为非cf投诉，我发现我无法分配x和y坐标。通过检查横截面示例中使用的netCDF文件，可以清楚地看到x和y数组表示具有GeoX和GeoY坐标轴类型的空间网格。然而，我不清楚的是如何为我的文件创建一个类似的数组。

纬度和经度数组不允许尝试手动折叠它们(如选择模式)(这是我尝试创建的一个合适的坐标向量)

如果有人能帮助解决这个问题，我将不胜感激。谢谢你！下面是我正在运行的代码。

import sys
#  pygrib is installed to the "herbie" environment, so we add the path
#  to this so we can load those packages.  Pygrib is used to get the 
#  complete set of keys for a message so we can try to see what is
#  available.
sys.path.append('/miniconda3/envs/herbie/lib/python3.8/site-packages')
import pygrib
%matplotlib inline
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.interpolate import cross_section
# NAM 218 AWIPS Grid - CONUS 
# (12-km Resolution; full complement of pressure level fields and some surface-based fields)
# https://nomads.ncep.noaa.gov/pub/data/nccf/com/nam/prod/nam.20210925/
filename = '/NAM/nam.t00z.awphys00.tm00.grib2'
#   NCEP GRIB files have non-unique keys, so pick from recommended list to filter
#    --Temperature and Pressure levels
#   Non-CF Compliant so don't use .parse_cf()
data = xr.open_dataset(filename, engine="cfgrib",filter_by_keys={'typeOfLevel': 'isobaricInhPa','name':'Temperature'})
#  We are missing a mapping projection, so add it:
data1=data.metpy.assign_crs(
grid_mapping_name='lambert_conformal_conic',
latitude_of_projection_origin=35,
longitude_of_central_meridian=-100,
standard_parallel=(30,60),
earth_radius=6371229.0)
print(data1)
#--------------------------  Print Output  ---------------------------------
#<xarray.Dataset>
#Dimensions:        (isobaricInhPa: 39, y: 428, x: 614)
#Coordinates:
#    time           datetime64[ns] 2021-09-25
#    step           timedelta64[ns] 00:00:00
#  * isobaricInhPa  (isobaricInhPa) float64 1e+03 975.0 950.0 ... 75.0 50.0
#    latitude       (y, x) float64 ...
#    longitude      (y, x) float64 ...
#    valid_time     datetime64[ns] 2021-09-25
#    metpy_crs      object Projection: lambert_conformal_conic
#Dimensions without coordinates: y, x
#Data variables:
#    t              (isobaricInhPa, y, x) float32 ...
#Attributes:
#    GRIB_edition:            2
#    GRIB_centre:             kwbc
#    GRIB_centreDescription:  US National Weather Service - NCEP 
#    GRIB_subCentre:          0
#    Conventions:             CF-1.7
#    institution:             US National Weather Service - NCEP 
#    history:                 2021-09-28T22:45 GRIB to CDM+CF via cfgrib-0.9.9...
#  Missing x- and y-coordinates
#----------------------------------------------------------------------
#  Try to add coordinates for y and x
data2=data1.metpy.assign_y_x()
#  Output Error:
# ValueError: Projected y and x coordinates cannot be collapsed to 1D within tolerance. Verify that your latitude and longitude coordinates correspond to your CRS coordinate.'''

Answer 1

所以这里的问题是由您包含的错误消息指出的:

验证经纬度坐标与CRS坐标对应

使用提供的CRS参数，计算出的x/y坐标似乎不是1D值，这使得它们在使用上存在问题，并表明投影参数通常不太正确。

我不确定您在assign_crs中使用的参数来自何处，但它们对于NAM 218网格不正确。您需要的参数是:

grid_mapping_name = "lambert_conformal_conic",
latitude_of_projection_origin = 25.0,
longitude_of_central_meridian = 265.0,
standard_parallel = 25.0,
earth_radius = 6371229.0

获得这些最可靠的方法是使用pygrib读取消息。假设一个文件中的所有GRIB消息都来自相同的坐标系(对于那些NWS文件来说很可能就是这种情况)，您可以这样做:

import pygrib
from pyproj import CRS
grib = pygrib.open(filename)
msg = grib.message(1)
data1 = data.metpy.assign_crs(CRS(msg.projparams).to_cf())
data2 = data1.metpy.assign_y_x()

来自pygrib的消息提供了适当的pyproj坐标，可以很容易地映射到必要的Cf参数。

如果cfgrib可以本地提供这些参数，那就太好了，但目前情况并非如此。这是一个开放的问题。