使用cartopy可视化Himawari-8 NetCDF数据的问题

Question

我正在尝试使用 Python 绘制 Himawari-8 AHI 数据（从 THREDD 下载）并遵循这个简单笔记本（转到 Himawari-8 部分）。第一次通过 imshow 绘制时一切顺利。

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs

# Data paths
h8b1_path = '0800_20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc'
h8b2_path = '0800_20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc'
h8b3_path = '0800_20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc'

# Open data
h8b1 = Dataset(h8b1_path)
h8b2 = Dataset(h8b2_path)
h8b3 = Dataset(h8b3_path)

for item in h8b1.dimensions:
    print("dimensions: {}, name: {}\n".format(h8b1.dimensions[item].name, 
                                              h8b1.dimensions[item].size))

vars = h8b1.variables.keys()
for item in vars:
    print('Variable:\t{}'.format(item))
    print('Dimensions:\t{}'.format(h8b1[item].dimensions))
    print('Shape:\t{}\n'.format(h8b1[item].shape))

b1 = h8b1.variables['channel_0001_scaled_radiance'][0,:,:]
b2 = h8b2.variables['channel_0002_scaled_radiance'][0,:,:]
b3 = h8b3.variables['channel_0003_scaled_radiance'][0,:,:]

x = h8b1.variables['x'][:]
y = h8b1.variables['y'][:]

b1 = h8b1.variables['channel_0001_scaled_radiance'][0,:,:]
b2 = h8b2.variables['channel_0002_scaled_radiance'][0,:,:]
b3 = h8b3.variables['channel_0003_scaled_radiance'][0,:,:]

x = h8b1.variables['x'][:]
y = h8b1.variables['y'][:]

fig = plt.figure(figsize=(12, 12))
spec = fig.add_gridspec(2, 4, hspace=0, wspace=0.1)

# Doing something here as to why I used gridspec.

ax0 = fig.add_subplot(spec[0, :])
ax0.imshow(rgb)

plt.show()

但是，当使用/应用来自 Cartopy 的 CRS 并重新绘制“Himawari 子集”时，它只显示一个空白图。

plt.figure(figsize=(12,12))

# Or perhaps the projection scheme being used???

# Define the projection information of the image
img_proj = ccrs.Geostationary(central_longitude=h8b1['geostationary'].longitude_of_projection_origin, 
                                  satellite_height=h8b1['geostationary'].satellite_height)

# The extents of the image we are plotting
img_extent = (X[0], X[-1], Y[-1], Y[0])

# Setup the axes projection
ax = plt.axes(projection=ccrs.Miller(central_longitude=float(h8b1['geostationary'].longitude_of_projection_origin)))

# Extent of the axes in the plot
ax.set_extent([117, 125, 13, 20]) # Area of Luzon Landmass

# Plot image
ax.imshow(rgb, transform=img_proj, extent=img_extent, origin='upper')

有办法解决这个问题吗？此处附有示例数据以及所附代码。

Answer 1

您没有展示如何创建您尝试绘制的

rgb

变量，这使得很难知道发生了什么。请记住，在使用 RGB(A) 定义时，matplotlib 对于值的范围是灵活的，但这种灵活性也可能会让您感到困扰。整数假定在 0-255 范围内，浮点数在 0-1 范围内。因此，如果您的 RGB（数组？）在 0-255 之间浮动，它可能会全部剪辑为白色。

肯定会出现问题的是单位不匹配。数据投影假定以米为单位，但 NetCDF 数据提供的 x/y 单位为公里。这将导致数据在投影原点附近绘制为小图像，可能也在您设置的地图范围之外。

下面的示例给出了解决此问题的示例。我改变了一些在某种程度上主观的东西。

• 使用 xarray 使读取数据稍微容易一些，但如果您想避免依赖性，NetCDF4 也应该可以正常工作。
• 在

  pcolorfast

  上使用

  imshow

  进行绘图。我倾向于避免使用 Cartopy 的

  imshow

  （不同数据与地图投影）对图像进行动态加热，尤其是像对地静止这样的投影，其中包含的点甚至不在地球上。这也是为什么我只读取一个子集（粗略地使用

  .isel

  ），这会加快速度，但也避免了重新投影时的那些“奇点”。

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pathlib import Path
import numpy as np

# I/O
data_dir = Path(r"C:\Temp")

with xr.open_mfdataset((
    data_dir / '0800_20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc',
    data_dir / '0800_20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc',
    data_dir / '0800_20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc',
)) as ds:

    ds = ds.isel(x=slice(1000, 4000), y=slice(1000, 3000), time=0)

    central_longitude = ds['geostationary'].longitude_of_projection_origin
    satellite_height = ds['geostationary'].satellite_height

    # convert km to m
    mapx = ds['x'].to_numpy() * 1000.
    mapy = ds['y'].to_numpy() * 1000.

    rgb = np.stack((
        ds['channel_0001_scaled_radiance'].to_numpy(),
        ds['channel_0002_scaled_radiance'].to_numpy(),
        ds['channel_0003_scaled_radiance'].to_numpy(),
    ), axis=-1)


# stretch RGB values a little
rgb_stretched = np.clip((rgb/0.7)**0.7, 0, 1)

# define projections
data_proj = ccrs.Geostationary(
    central_longitude=central_longitude, 
    satellite_height=satellite_height,
)
map_proj =  ccrs.Miller(central_longitude=central_longitude)

# plotting
fig, ax = plt.subplots(
    figsize=(6, 6), facecolor="w", dpi=100, layout="compressed",
    subplot_kw=dict(projection=map_proj),
)

pcm = ax.pcolorfast(mapx, mapy, rgb_stretched, transform=data_proj)

ax.coastlines(color="w")
ax.set_extent([115, 127, 11, 22])

为了可视化目的，我稍微缩放了原始 RGB，并且还设置了比您稍宽的范围以显示更多上下文。

对于您拥有的相对较小的子集，我可能不会将地图投影更改为与数据投影不同的东西。这只是引入了额外的重采样，而没有太多额外的好处（对我来说），当然这样做是完全有效的。