绘制 Himawari-8 NetCDF 数据 - 轴与数组形状不匹配

Question

这和去年的这篇post很有关系，可以参考。我正在尝试可视化以菲律宾为中心的 Himawari-8 AHI，并通过 NCI Thredds Server 下载。不幸的是，我遇到了一些问题，特别是在制作真彩色图像时（类似于这个），即通过 Cartopy 堆叠和绘制它。

使用 xarray 读取时，波段 1、2 和 4（素食）具有相似的 1000 m 分辨率和尺寸，但波段 3 的分辨率为 500 m。鉴于它们具有不同的分辨率，我最初对每个 NetCDF 数据的 x 和 y 坐标进行了切片，然后将它们与

xr.merge

合并。

我尝试用

np.stack

堆叠每个带，但不幸的是，它给出了一个错误；

Axes don't match array shape. Got (4, 1), expected (2999, 2999).

目前代码运行如下：

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("New folder (3)")

ds1 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc')
ds2 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc')
ds3 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_500-HIMAWARI8-AHI.nc')
ds4 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B04-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc')

# x and y coordinates are in meters already. 
# We select a portion of the FDK focused around Luzon i.e. Central to Northern Luzon

dx1 = ds1.isel(x=slice(3000, 4000), y=slice(3000, 4000))
#dx1

dx2 = ds2.isel(x=slice(3000, 4000), y=slice(3000, 4000))
#dx2

# Has a different resolution, matched up to the 1000 m of other bands
dx3 = ds3.isel(x=slice(6000, 8000), y=slice(6000, 8000))
#dx3

dx4 = ds4.isel(x=slice(3000, 4000), y=slice(3000, 4000))
#dx4

dx = xr.merge([dx1, dx2, dx3, dx4])

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

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

rgb = np.stack((
    dx['channel_0001_scaled_radiance'].to_numpy(), # Band 1 is blue (0.47 um)
    dx['channel_0002_scaled_radiance'].to_numpy(), # Band 2 is green (0.51 um)
    dx['channel_0003_scaled_radiance'].to_numpy(), # Band 3 is red (0.64 um)
    dx['channel_0004_scaled_radiance'].to_numpy(), # Band 4 is "veggie" infrared (0.865 um)
))

# stretch RGB values a little
rgb_stretched = np.clip((rgb/0.85)**0.85, 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=(20, 12), facecolor="w", dpi=300,
    subplot_kw=dict(projection=map_proj),
)

datacrs = ccrs.PlateCarree()

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

ax.coastlines(color="w")

我相信由于我的Python新手技能，我在这个过程中错过了或者做了一些不正确的事情。有没有办法解决这个问题，以便可以显示 Himawari-8 卫星图像？这里附上代码中使用的示例数据。

Answer 1

这有几个问题。但首先我要提一下，有专门的库可用于此目的，例如 SatPy (https://satpy.readthedocs.io/en/stable/index.html#)，这些对某些人来说可能更方便。

问题 1）与 Xarray 的合并是在（数据）坐标上执行的，对于 500m 波段来说，这些坐标是不同的（偏移）。这引入了很多“介于两者之间”的 NaN 值。

解决方案 1) 您可以从高分辨率带中选择每隔一个像素，假设两个网格都是一条线（例如

.slice(6000, 8000, 2)

）。但更稳健的方法是将所有数据插值/重新采样到一个公共网格。最简单的方法是将频段 3 重新采样到其他频段，但您也可以以相反的方式进行，以保留频段 3 的 500m 分辨率。所以：

dx = xr.merge([dx1, dx2, dx3.interp_like(dx1), dx4])

问题 2) Matplotlib 需要 RGB 格式，并在最后一个维度（x、y、颜色）上包含颜色信息。您的示例还具有时间维度，您有（颜色，时间，x，y）。

解决方案 2）由于您只有一个时间戳，因此您可以通过将

isel(time=0)

添加到您的初始切片来删除它。这降低了时间维度。例如对于第一个乐队：

dx1 = ds1.isel(time=0, x=slice(3000, 4000, 1), y=slice(3000, 4000, 1))

并且在堆叠 Numpy 数组时指定最后一个维度会将颜色信息添加到最后一个维度（而不是默认的第一个维度），因此：

rgb = np.stack((...), axis=-1)

我还会删除第四个 NIR 波段，否则这将被解释为 alpha（透明度），我怀疑这不是您想要的，但请随意将其添加回来，或者为最终图像选择不同的波段组合。

问题 3）Matplotlib 的

pcolorfast

（或

pcolormesh

等）期望坐标数组来指定边缘，因此与数据（N+1，M+1）相比，它需要一个额外的坐标元素。这可能是 Matplotlib 最近的变化（？）。

解决方案 3）由于您正在处理常规网格，因此更快的解决方案是仅提供网格的外边缘/范围。这可以通过选择定义该像素中心的外部坐标，并添加/减去一半的分辨率（在下面的示例中为 500m）来完成。

mapx = (mapx[0]-500, mapx[-1]+500)
mapy = (mapy[0]+500, mapy[-1]-500)

然后给出：

完整的代码，为了完整性。可能可以将 Xarray 部分重写得更干净一些，但这仍然接近 OP 的原始代码。

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

data_dir = Path("D:\Temp\him8")

ds1 = xr.open_dataset(data_dir / '20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc')
ds2 = xr.open_dataset(data_dir / '20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc')
ds3 = xr.open_dataset(data_dir / '20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_500-HIMAWARI8-AHI.nc')

dx1 = ds1.isel(time=0, x=slice(3000, 4000, 1), y=slice(3000, 4000, 1))
dx2 = ds2.isel(time=0, x=slice(3000, 4000, 1), y=slice(3000, 4000, 1))
dx3 = ds3.isel(time=0, x=slice(6000, 8000, 1), y=slice(6000, 8000, 1))

dx = xr.merge([dx1, dx2, dx3.interp_like(dx1)])

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

mapx = dx1['x'].to_numpy()
mapy = dx1['y'].to_numpy()

mapx = (mapx[0]-500, mapx[-1]+500)
mapy = (mapy[0]+500, mapy[-1]-500)

rgb = np.stack((
    dx['channel_0003_scaled_radiance'].to_numpy(), # Band 3 is red (0.64 um)
    dx['channel_0002_scaled_radiance'].to_numpy(), # Band 2 is green (0.51 um
    dx['channel_0001_scaled_radiance'].to_numpy(), # Band 1 is blue (0.47 um)
), axis=-1)

rgb_stretched = np.clip((rgb/0.85)**0.85, 0, 1)

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

fig, ax = plt.subplots(
    figsize=(10, 10), facecolor="w", dpi=100,
    subplot_kw=dict(projection=map_proj),
)

pcm = ax.pcolorfast(mapx, mapy, rgb_stretched, transform=data_proj)
ax.coastlines(color="w")
ax.set_extent((*mapx, *mapy), crs=data_proj)