将 GOES-16 地球静止数据投影到 Plate Carree Cartopy

Projecting GOES-16 Geostationary data into Plate Carree Cartopy

我拼命尝试将一些地球静止数据从 GOES-16 netCDF 文件投影到不同的投影。我可以获取要重新投影的背景图,但似乎无法获取要遵循的数据。

我还不是很精通这方面,但这是我目前所掌握的:

正在通过NetCDF4读取数据:

from netCDF4 import Dataset
nc = Dataset('OR_ABI-L1b-RadF- 
    M3C13_G16_s20182831030383_e20182831041161_c20182831041217.nc')

data = nc.variables['Rad'][:]

我在这里尝试获取地球同步信息:

sat_h = nc.variables['goes_imager_projection'].perspective_point_height
X = nc.variables['x'][:] * sat_h
Y = nc.variables['y'][:] * sat_h

# Satellite longitude
sat_lon = 
    nc.variables['goes_imager_projection'].longitude_of_projection_origin

# Satellite sweep
sat_sweep = nc.variables['goes_imager_projection'].sweep_angle_axis

这里我从 .nc 文件中获取投影数据:

proj_var = nc.variables['goes_imager_projection']

sat_height = proj_var.perspective_point_height
central_lon = proj_var.longitude_of_projection_origin
semi_major = proj_var.semi_major_axis
semi_minor = proj_var.semi_minor_axis
print proj_var

<type 'netCDF4._netCDF4.Variable'>
int32 goes_imager_projection()
    long_name: GOES-R ABI fixed grid projection
    grid_mapping_name: geostationary
    perspective_point_height: 35786023.0
    semi_major_axis: 6378137.0
    semi_minor_axis: 6356752.31414
    inverse_flattening: 298.2572221
    latitude_of_projection_origin: 0.0
    longitude_of_projection_origin: -75.0
    sweep_angle_axis: x
unlimited dimensions: 
current shape = ()
filling on, default _FillValue of -2147483647 used

这是我的一小段相关代码:

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

globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)
proj = ccrs.Geostationary(central_longitude=central_lon, 
    satellite_height=sat_height, globe=globe)

ax = fig.add_subplot(1, 1, 1, projection=proj)

IR_img = ax.imshow(data[:,:],origin='upper',extent=(X.min(), X.max(), Y.min(), Y.max()),
    cmap=IR_cmap,interpolation='nearest',vmin=162.,vmax=330.)

还有一张大家打得很好的照片: Data and map working

当我尝试说 Plate Carree 投影时,我尝试:

proj = ccrs.PlateCarree(central_longitude=central_lon,globe=globe)

还有一张我失败的照片: Data and map not working

我试过在 imshow 方法中弄乱范围,我试过添加

transform=proj 

在 imshow 中,没有运气,它只是挂起,我必须重新启动内核。

很明显是我理解不够。如果有人可以快速轻松地 help/explain 我想改变地球同步投影的方式,我将不胜感激。

我 运行 过时 python2。

感谢观看。

编辑:多亏了 DopplerShift 和 ajdawson 的洞察力,问题似乎得到了解决,我想我可能有点 impatient/ignorant 全盘转换需要多长时间。

看来您需要指定 transform 关键字才能显示。这个关键字告诉 cartopy 你的数据在什么坐标中,在这种情况下应该是地球静止的。

我没有你的数据集,所以我无法对此进行测试,但下面的代码片段说明了这个概念。投影和变换是独立的,因此您应该同时定义两者。 transform 参数的值(下例中的 crs)对于数据集是固定的,但投影可以是任何你喜欢的(包括与 crs 相同)。

查看这个重新投影地球静止图像的示例:https://scitools.org.uk/cartopy/docs/v0.16/gallery/geostationary.html#sphx-glr-gallery-geostationary-py. Also see the guide to projection and transform arguments here: https://scitools.org.uk/cartopy/docs/v0.16/tutorials/understanding_transform.html

globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)
crs = ccrs.Geostationary(central_longitude=central_lon, 
                         satellite_height=sat_height, globe=globe)
proj = ccrs.PlateCarree(central_longitude=central_lon, globe=globe)

ax = fig.add_subplot(1, 1, 1, projection=proj)

IR_img = ax.imshow(data[:,:], origin='upper',
                   extent=(X.min(), X.max(), Y.min(), Y.max()),
                   transform=crs,
                   cmap=IR_cmap,
                   interpolation='nearest', vmin=162., vmax=330.)