来自 Grib2 文件的 Cartopy 图像与同一 CRS 中的海岸线和边界不对齐

Question

这看起来会很简单，可以修复我的代码，但我认为我目前对代码看得太多了，需要重新审视它。我只是想引入我从 NCEP 下载的用于 HRRR 模型的 Grib2 文件。根据他们的信息，网格类型是 Lambert Conformal，角的纬度范围为 (21.13812, 21.14055, 47.84219, 47.83862)，角的经度范围为 (-122.7195, -72.28972, -60.91719, -134.0955)模型域。

在尝试放大我感兴趣的区域之前，我只想简单地在适当的 CRS 中显示图像，但是当我尝试为模型的域执行此操作时，我得到的边界和海岸线落入其中范围，但从 Grib2 文件生成的实际图像只是放大了。我尝试使用 extent=[my domain extent] 但它似乎总是让我正在测试它的笔记本崩溃。这是我的代码和我从中获得的相关图像它。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
from mpl_toolkits.basemap import Basemap
from osgeo import gdal
gdal.SetConfigOption('GRIB_NORMALIZE_UNITS', 'NO')

plt.figure()
filename='C:\Users\Public\Documents\GRIB\hrrr.t18z.wrfsfcf00.grib2'
grib = gdal.Open(filename, gdal.GA_ReadOnly)

z00 = grib.GetRasterBand(47)
meta00 = z00.GetMetadata()
band_description = z00.GetDescription()

bz00 = z00.ReadAsArray()

latitude_south = 21.13812 #38.5
latitude_north = 47.84219 #50
longitude_west = -134.0955 #-91
longitude_east = -60.91719 #-69

fig = plt.figure(figsize=(20, 20))

title= meta00['GRIB_COMMENT']+' at '+meta00['GRIB_SHORT_NAME']
fig.set_facecolor('white')

ax = plt.axes(projection=ccrs.LambertConformal())

ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.coastlines(resolution='110m')

ax.imshow(bz00,origin='upper',transform=ccrs.LambertConformal())

plt.title(title)
plt.show()

Returns Just Grib File

如果我改变：

ax = plt.axes(projection=ccrs.LambertConformal()

到

ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-95.5, 
central_latitude=38.5,cutoff=21.13)

我得到了我的边界，但我的实际数据没有对齐，它创建了我正在配音的蝙蝠侠情节。

Batman Plot

即使我放大域并且仍然有我的边框，也会出现类似的问题。 Grib 文件中的基础数据不会更改以对应于我要获取的内容。

因此，正如我已经说过的，这可能是我所缺少的一个简单的修复方法，但如果不是，那么很高兴知道我在搞砸哪个步骤或哪个过程，我可以吸取教训，以后不做！

更新 1： 我已经添加并更改了一些代码，然后又回到了只显示图像而不显示边界和海岸线的状态。

test_extent = [longitude_west,longitude_east,latitude_south,latitude_north]
ax.imshow(bz00,origin='upper',extent=test_extent)

这给了我下面的图像。 Looks exactly like image 1.

我注意到的另一件事可能是所有这一切的根本原因是，当我打印出 plt.gca().get_ylim() 和 plt.gca().get_xlim() 的值时，我变得非常不同值取决于显示的内容。

Answer 1

我的问题似乎是因为 Grib 文件无论是否可以在其他程序中正确显示，都不能很好地与开箱即用的 Matplotlib 和 Cartopy 一起播放。或者至少与我使用的 Grib 文件无关。为此，您可以从 NCEP HRRR 模型中获得 here or here。

如果将文件从 Grib2 格式转换为 NetCDF 格式，一切似乎都能正常工作，并且我能够在地图上获得我想要的边界、海岸线等内容。我附上了下面的代码和输出以显示它是如何工作的。此外，我还手动选择了一个我想显示的数据集，以与我以前的代码进行测试，因此如果您想查看文件中可用的其余数据集，您需要使用 ncdump 或类似的东西来查看数据集上的信息.

import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
import cartopy.feature as cfeature
from osgeo import gdal
gdal.SetConfigOption('GRIB_NORMALIZE_UNITS', 'NO')

nc_f = 'C:\Users\Public\Documents\GRIB\test.nc'  # Your filename
nc_fid = Dataset(nc_f, 'r')  # Dataset is the class behavior to open the 
                             # file and create an instance of the ncCDF4 
                             # class

# Extract data from NetCDF file
lats = nc_fid.variables['gridlat_0'][:]  
lons = nc_fid.variables['gridlon_0'][:]
temp = nc_fid.variables['TMP_P0_L1_GLC0'][:]

fig = plt.figure(figsize=(20, 20))
states_provinces = cfeature.NaturalEarthFeature(category='cultural', \
      name='admin_1_states_provinces_lines',scale='50m', facecolor='none')
proj = ccrs.LambertConformal()
ax = plt.axes(projection=proj)

plt.pcolormesh(lons, lats, temp, transform=ccrs.PlateCarree(), 
              cmap='RdYlBu_r', zorder=1)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':', zorder=2)
ax.add_feature(states_provinces, edgecolor='black')
ax.coastlines()

plt.show()

地图最终预览

来自 Grib2 文件的 Cartopy 图像与同一 CRS 中的海岸线和边界不对齐

Cartopy Image from Grib2 file does not align with Coastlines and Borders in same CRS

python

matplotlib

grib

cartopy