绘制Himawari-8 NetCDF数据时轴与数组形状不匹配的问题求助

我参考去年的相关帖子，尝试用Cartopy对NCI Thredds Server下载的菲律宾中心Himawari-8 AHI数据制作真彩色图像。

用xarray读取数据后发现：Band1、2、4（植被波段）分辨率为1000m且维度一致，但Band3分辨率为500m。针对分辨率差异，我先对各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卫星图像。

内容的提问来源于stack exchange，提问作者Senkei