这似乎是一个简单的事情,可以修复我的代码,但我想我现在已经对代码看了太多了,需要重新审视它。我只是尝试导入从 NCEP 下载的 HRRR 模型的 Grib2 文件。根据他们的信息,网格类型为兰伯特等角,角点纬度范围为(21.13812, 21.14055, 47.84219, 47.83862),角点经度范围为(-122.7195, -72.28972, -60.91719, -134.0955)模型域。
在尝试放大我感兴趣的区域之前,我只想简单地在适当的 CRS 中显示图像,但是当我尝试对模型的域执行此操作时,我得到的边界和海岸线都在该范围内,但是从 Grib2 文件生成的实际图像只是被放大。我尝试使用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()
如果我改变:
ax = plt.axes(projection=ccrs.LambertConformal()
至
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-95.5,
central_latitude=38.5,cutoff=21.13)
我得到了边框,但我的实际数据没有对齐,它创建了我所说的 bat 侠情节。
即使我放大域并且边框仍然存在,也会出现类似的问题。 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()
根据显示的内容,我得到的值有很大不同。
最佳答案
看来我的问题源于这样一个事实:Grib 文件无论是否可以在其他程序中正确显示,都不能很好地与 Matplotlib 和 Cartopy 一起使用。或者至少我使用的 Grib 文件没有。为此,您可以从 NCEP HRRR 模型中获取 here或here .
如果将文件从 Grib2 格式转换为 NetCDF 格式,一切似乎都会很好地工作,并且我能够在 map 上获得我想要的边界、海岸线等内容。我已附上代码和下面的输出来展示它是如何工作的。此外,我还手动选择了一个要显示的数据集,以与之前的代码进行比较,因此如果您想查看文件中可用的其余数据集,则需要使用 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()
map 最终预览
关于python - Grib2 文件中的 Cartopy 图像与同一 CRS 中的 CoaSTLines 和边框不对齐,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/48393305/