使用 Metpy assign_y_x 坐标和手动 CRS 时坐标对齐关闭
Coordinate alignment off when using Metpy assign_y_x coordinates with manual CRS
我在使用为网格映射分配 y 和 x 坐标值的函数时遇到问题。本质上，我的最终目标是能够从一组给定的 lat/lon 对中提取网格值，这似乎比它应该的要复杂得多（也许我的方法不是最好的）。我正在使用标准的 RTMA grib 文件并使用 xarray 和 cfgrib 引擎打开。 RTMA 数据没有附加到数据的网格映射，因此需要手动添加 CRS。我试图通过直接从 grib 文件获取信息并使用 pygrib 进行确认来提供它。问题是 x/y 坐标值似乎没有与它们各自的 lat/lon 坐标正确对齐，因此当尝试为点提取建立索引时，您可以关闭一些网格单元格（可以有与降水梯度有很大差异）。以下是我到目前为止所做的工作：
rtma = xr.open_dataset("/rtma2p5.2022022023z.conus.grb2", engine="cfgrib")
print(rtma)
<xarray.Dataset>
Dimensions: (y: 1377, x: 2145)
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
heightAboveGround float64 ...
latitude (y, x) float64 ...
longitude (y, x) float64 ...
valid_time datetime64[ns] ...
surface float64 ...
Dimensions without coordinates: y, x
Data variables:
t2m (y, x) float32 ...
d2m (y, x) float32 ...
sp (y, x) float32 ...
vis (y, x) float32 ...
Attributes:
GRIB_edition: 2
GRIB_centre: kwbc
GRIB_centreDescription: US National Weather Service - NCEP
GRIB_subCentre: 4
Conventions: CF-1.7
institution: US National Weather Service - NCEP
history: 2022-04-19T15:21 GRIB to CDM+CF via cfgrib-0.9.1...
我们有纬度和经度坐标，没有购买 x/y 坐标值，即使在尝试使用 metpy.parse_cf() 之后也是如此。所以，是时候手动添加网格映射了。
使用 pygrib 我们可以确认一些 CRS 值
grib = pygrib.open("/rtma2p5.2022022023z.conus.grb2")
msg = grib.message(1)
CRS(msg.projparams).to_cf()
{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371200,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",25,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",265,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",25,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",25,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
'semi_major_axis': 6371200.0,
'semi_minor_axis': 6371200.0,
'inverse_flattening': 0.0,
'reference_ellipsoid_name': 'unknown',
'longitude_of_prime_meridian': 0.0,
'prime_meridian_name': 'Greenwich',
'geographic_crs_name': 'unknown',
'horizontal_datum_name': 'unknown',
'projected_crs_name': 'unknown',
'grid_mapping_name': 'lambert_conformal_conic',
'standard_parallel': (25.0, 25.0),
'latitude_of_projection_origin': 25.0,
'longitude_of_central_meridian': 265.0,
'false_easting': 0.0,
'false_northing': 0.0}
很好，现在我们将分配网格映射以获得 x/y 坐标值。我们将使用 grib 文件中的属性数据来添加它，因为我们使用 pygrib
确认它看起来正确
# Assign the grid mapping
grid = rt.metpy.assign_crs({
"semi_major_axis": 6371200.0,
"semi_minor_axis": 6371200.0,
"grid_mapping_name": "lambert_conformal_conic",
"standard_parallel": [
rt["t2m"].attrs["GRIB_Latin1InDegrees"],
rt["t2m"].attrs["GRIB_Latin2InDegrees"]
],
"latitude_of_projection_origin": rt["t2m"].attrs["GRIB_LaDInDegrees"],
"longitude_of_central_meridian": rt["t2m"].attrs["GRIB_LoVInDegrees"],
}).metpy.assign_y_x()
print(grid)
<xarray.Dataset>
Dimensions: (y: 1377, x: 2145)
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
heightAboveGround float64 ...
latitude (y, x) float64 20.19 20.2 20.2 ... 50.12 50.11 50.11
longitude (y, x) float64 238.4 238.5 238.5 ... 299.1 299.1 299.1
valid_time datetime64[ns] ...
surface float64 ...
metpy_crs object Projection: lambert_conformal_conic
* y (y) float64 -2.638e+05 -2.612e+05 ... 3.228e+06 3.231e+06
* x (x) float64 -2.763e+06 -2.761e+06 ... 2.679e+06 2.682e+06
Data variables:
t2m (y, x) float32 ...
d2m (y, x) float32 ...
sp (y, x) float32 ...
vis (y, x) float32 ...
Attributes:
GRIB_edition: 2
GRIB_centre: kwbc
GRIB_centreDescription: US National Weather Service - NCEP
GRIB_subCentre: 4
Conventions: CF-1.7
institution: US National Weather Service - NCEP
history: 2022-04-19T15:21 GRIB to CDM+CF via cfgrib-0.9.1...
当尝试从 lat/lon 对中提取值时，问题就来了。我可以将示例 lat/lon 对转换为网格投影中的坐标：
grid_crs = ccrs.LambertConformal(central_longitude = grid["t2m"].attrs["GRIB_LoVInDegrees"],
central_latitude = grid["t2m"].attrs["GRIB_LaDInDegrees"],
standard_parallels = [
grid["t2m"].attrs["GRIB_Latin1InDegrees"],
grid["t2m"].attrs["GRIB_Latin2InDegrees"]
]
)
lat = 35.36234
lon = 283.28369
x_t, y_t = grid_crs.transform_point(lon, lat, src_crs=ccrs.PlateCarree())
然后使用
grid.sel(x=x_t, y=y_t, method='nearest')
或
grid.interp(x=x_t, y=y_t, method='nearest')
但是我注意到这些值有点偏差。经过大量的试验和错误并绘制完成插值或最近邻点的值后，我注意到这些点已经关闭。所以我回去绘制边缘值（grib 文件中的第一个 x/y 坐标值对）并注意到它没有完全对齐。
print(grid.coords["x"][0].values)
-2763204.4992319937
print(grid.coords["y"][0].values)
-263789.4687076026
绘制网格的pcolormesh以及网格中的lat/lon值后，可以看到第一个x/y坐标点（红色）与左下角不匹配（蓝色）lat/lon 对。当尝试索引或插值点值时，这会导致整个网格出现错误值。
最后，这张图片包含很多内容，但希望它能帮助您理解。
例如，此处绘制的是 t2m 字段的 pcolormesh。浅蓝色点是文件中的 lat/lon 坐标。黑环是我尝试从中提取数据的请求位置。紫色的点是车站周围 lat/lon 坐标的 3x3 网格（我使用此处描述的方法找到了离车站最近的索引：）。橙色点是相同的 3x3 网格，但使用 x/y 坐标值对绘制。这表明 x/y 坐标值对在绘制时不对应于相同的纬度和经度。发生这种情况时，您 运行 进入问题索引点。当您使用 nearest 方法 运行 .sel 时，您会看到它抓取了最近的 x/y 坐标对（亮绿色点与橙色点重叠），它对应于更远的 lat/lon 对离开（亮绿色实心点）。尝试插值时可以看到同样的情况；它可以正确定位车站（粉红点与黑点重叠），但它对应的 lat/lon 值又离得更远，给你一个不同的值。
希望这是有道理的，它有点难以解释，但希望最后一张图片有所帮助。有什么办法可以解决这个问题吗？或者也许是一种更好的提取点值的方法？
造成这种转变的最可能原因是数据集中使用的 CRS（对于给定的 latitude/longitudes 和计算的 y/x）与用于从所需的 latitude/longitude 对计算投影 y/x 对。特别要注意的是，数据集的 CRS 是用以下内容构建的：
"semi_major_axis": 6371200.0,
"semi_minor_axis": 6371200.0,
但是在使用 ccrs.LambertConformal() 构造 grid_crs 时没有提供椭球参数，这意味着使用了默认的椭球 (wgs84)。
幸运的是，MetPy 具有处理数据集中 CRS 的实用程序，应该可以更轻松地保持 CRS 和全球一致。我建议使用 MetPy 的访问器，而不是自己构建网格 CRS，并保留 lat/lon CRS 未指定：
grid_crs = grid["t2m"].metpy.cartopy_crs
latlon_crs = ccrs.PlateCarree(globe=grid["t2m"].metpy.cartopy_globe)
x_t, y_t = grid_crs.transform_point(lon, lat, src_crs=latlon_crs)
我在使用为网格映射分配 y 和 x 坐标值的函数时遇到问题。本质上，我的最终目标是能够从一组给定的 lat/lon 对中提取网格值，这似乎比它应该的要复杂得多（也许我的方法不是最好的）。我正在使用标准的 RTMA grib 文件并使用 xarray 和 cfgrib 引擎打开。 RTMA 数据没有附加到数据的网格映射，因此需要手动添加 CRS。我试图通过直接从 grib 文件获取信息并使用 pygrib 进行确认来提供它。问题是 x/y 坐标值似乎没有与它们各自的 lat/lon 坐标正确对齐，因此当尝试为点提取建立索引时，您可以关闭一些网格单元格（可以有与降水梯度有很大差异）。以下是我到目前为止所做的工作：
rtma = xr.open_dataset("/rtma2p5.2022022023z.conus.grb2", engine="cfgrib")
print(rtma)
<xarray.Dataset>
Dimensions: (y: 1377, x: 2145)
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
heightAboveGround float64 ...
latitude (y, x) float64 ...
longitude (y, x) float64 ...
valid_time datetime64[ns] ...
surface float64 ...
Dimensions without coordinates: y, x
Data variables:
t2m (y, x) float32 ...
d2m (y, x) float32 ...
sp (y, x) float32 ...
vis (y, x) float32 ...
Attributes:
GRIB_edition: 2
GRIB_centre: kwbc
GRIB_centreDescription: US National Weather Service - NCEP
GRIB_subCentre: 4
Conventions: CF-1.7
institution: US National Weather Service - NCEP
history: 2022-04-19T15:21 GRIB to CDM+CF via cfgrib-0.9.1...
我们有纬度和经度坐标，没有购买 x/y 坐标值，即使在尝试使用 metpy.parse_cf() 之后也是如此。所以，是时候手动添加网格映射了。
使用 pygrib 我们可以确认一些 CRS 值
grib = pygrib.open("/rtma2p5.2022022023z.conus.grb2")
msg = grib.message(1)
CRS(msg.projparams).to_cf()
{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371200,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",25,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",265,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",25,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",25,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
'semi_major_axis': 6371200.0,
'semi_minor_axis': 6371200.0,
'inverse_flattening': 0.0,
'reference_ellipsoid_name': 'unknown',
'longitude_of_prime_meridian': 0.0,
'prime_meridian_name': 'Greenwich',
'geographic_crs_name': 'unknown',
'horizontal_datum_name': 'unknown',
'projected_crs_name': 'unknown',
'grid_mapping_name': 'lambert_conformal_conic',
'standard_parallel': (25.0, 25.0),
'latitude_of_projection_origin': 25.0,
'longitude_of_central_meridian': 265.0,
'false_easting': 0.0,
'false_northing': 0.0}
很好，现在我们将分配网格映射以获得 x/y 坐标值。我们将使用 grib 文件中的属性数据来添加它，因为我们使用 pygrib
确认它看起来正确# Assign the grid mapping
grid = rt.metpy.assign_crs({
"semi_major_axis": 6371200.0,
"semi_minor_axis": 6371200.0,
"grid_mapping_name": "lambert_conformal_conic",
"standard_parallel": [
rt["t2m"].attrs["GRIB_Latin1InDegrees"],
rt["t2m"].attrs["GRIB_Latin2InDegrees"]
],
"latitude_of_projection_origin": rt["t2m"].attrs["GRIB_LaDInDegrees"],
"longitude_of_central_meridian": rt["t2m"].attrs["GRIB_LoVInDegrees"],
}).metpy.assign_y_x()
print(grid)
<xarray.Dataset>
Dimensions: (y: 1377, x: 2145)
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
heightAboveGround float64 ...
latitude (y, x) float64 20.19 20.2 20.2 ... 50.12 50.11 50.11
longitude (y, x) float64 238.4 238.5 238.5 ... 299.1 299.1 299.1
valid_time datetime64[ns] ...
surface float64 ...
metpy_crs object Projection: lambert_conformal_conic
* y (y) float64 -2.638e+05 -2.612e+05 ... 3.228e+06 3.231e+06
* x (x) float64 -2.763e+06 -2.761e+06 ... 2.679e+06 2.682e+06
Data variables:
t2m (y, x) float32 ...
d2m (y, x) float32 ...
sp (y, x) float32 ...
vis (y, x) float32 ...
Attributes:
GRIB_edition: 2
GRIB_centre: kwbc
GRIB_centreDescription: US National Weather Service - NCEP
GRIB_subCentre: 4
Conventions: CF-1.7
institution: US National Weather Service - NCEP
history: 2022-04-19T15:21 GRIB to CDM+CF via cfgrib-0.9.1...
当尝试从 lat/lon 对中提取值时，问题就来了。我可以将示例 lat/lon 对转换为网格投影中的坐标：
grid_crs = ccrs.LambertConformal(central_longitude = grid["t2m"].attrs["GRIB_LoVInDegrees"],
central_latitude = grid["t2m"].attrs["GRIB_LaDInDegrees"],
standard_parallels = [
grid["t2m"].attrs["GRIB_Latin1InDegrees"],
grid["t2m"].attrs["GRIB_Latin2InDegrees"]
]
)
lat = 35.36234
lon = 283.28369
x_t, y_t = grid_crs.transform_point(lon, lat, src_crs=ccrs.PlateCarree())
然后使用
grid.sel(x=x_t, y=y_t, method='nearest')
或
grid.interp(x=x_t, y=y_t, method='nearest')
但是我注意到这些值有点偏差。经过大量的试验和错误并绘制完成插值或最近邻点的值后，我注意到这些点已经关闭。所以我回去绘制边缘值（grib 文件中的第一个 x/y 坐标值对）并注意到它没有完全对齐。
print(grid.coords["x"][0].values)
-2763204.4992319937
print(grid.coords["y"][0].values)
-263789.4687076026
绘制网格的pcolormesh以及网格中的lat/lon值后，可以看到第一个x/y坐标点（红色）与左下角不匹配（蓝色）lat/lon 对。当尝试索引或插值点值时，这会导致整个网格出现错误值。
最后，这张图片包含很多内容，但希望它能帮助您理解。
例如，此处绘制的是 t2m 字段的 pcolormesh。浅蓝色点是文件中的 lat/lon 坐标。黑环是我尝试从中提取数据的请求位置。紫色的点是车站周围 lat/lon 坐标的 3x3 网格（我使用此处描述的方法找到了离车站最近的索引：
希望这是有道理的，它有点难以解释，但希望最后一张图片有所帮助。有什么办法可以解决这个问题吗？或者也许是一种更好的提取点值的方法？
造成这种转变的最可能原因是数据集中使用的 CRS（对于给定的 latitude/longitudes 和计算的 y/x）与用于从所需的 latitude/longitude 对计算投影 y/x 对。特别要注意的是，数据集的 CRS 是用以下内容构建的：
"semi_major_axis": 6371200.0,
"semi_minor_axis": 6371200.0,
但是在使用 ccrs.LambertConformal() 构造 grid_crs 时没有提供椭球参数，这意味着使用了默认的椭球 (wgs84)。
幸运的是，MetPy 具有处理数据集中 CRS 的实用程序，应该可以更轻松地保持 CRS 和全球一致。我建议使用 MetPy 的访问器，而不是自己构建网格 CRS，并保留 lat/lon CRS 未指定：
grid_crs = grid["t2m"].metpy.cartopy_crs
latlon_crs = ccrs.PlateCarree(globe=grid["t2m"].metpy.cartopy_globe)
x_t, y_t = grid_crs.transform_point(lon, lat, src_crs=latlon_crs)