首页 文章

来自Grib2文件的Cartopy图像与同一CRS中的Coastlines和Borders不对齐

提问于
浏览
1

这看起来很简单,会修复我的代码,但我想我现在只是看了太多代码,需要对它进行一些新的研究 . 我只是想引入一个我从NCEP为HRRR模型下载的Grib2文件 . 根据他们的信息,网格类型是Lambert Conformal,角度的纬度为(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()

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文件的基础数据不会改变以符合我想要获得的内容 .

所以我已经说过这可能是一个很容易解决的问题,我只是缺少但是如果没有,那么知道我正在搞砸了哪个步骤或哪个过程我可以从中学习以便我不要在将来这样做!

Updated 1: 我已经添加并更改了一些代码,我回过头来只显示图像,而不显示边框和海岸线 .

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

这给了我以下图像 . Looks exactly like image 1.

另一件事是'm noticing which maybe the root cause of all of this is that when I' m打印出 plt.gca().get_ylim()plt.gca().get_xlim() 的值我根据显示的内容得到了非常不同的值 .

1 回答

  • 0

    似乎我的问题源于Grib文件,无论它是否可以在其他程序中正确显示,只是不能很好地与Matplotlib和Cartopy开箱即用 . 或者至少不使用我正在使用的Grib文件 . 为了这个可能帮助将来的其他人来自NCEP HRRR模型,你可以得到herehere .

    如果您将文件从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 的最终预览

    enter image description here

相关问题