基于xarray Dataset绘制卫星图像出现灰色方块，求正确方法

问题说明 #

现有xarray卫星数据集结构如下：

<xarray.Dataset>
Dimensions:                      (time: 1, y_geostationary: 4176,
                                  x_geostationary: 5548, variable: 1)
Coordinates:
  * time                         (time) datetime64[ns] 2021-01-01
  * variable                     (variable) <U3 'HRV'
  * x_geostationary              (x_geostationary) float64 -2.733e+06 ... 2.8...
  * y_geostationary              (y_geostationary) float64 1.395e+06 ... 5.57...
Data variables:
    data                         (time, y_geostationary, x_geostationary, variable) float16 dask.array<chunksize=(1, 100, 100, 1), meta=np.ndarray>
    x_geostationary_coordinates  (time, x_geostationary) float64 dask.array<chunksize=(1, 5548), meta=np.ndarray>
    y_geostationary_coordinates  (time, y_geostationary) float64 dask.array<chunksize=(1, 4176), meta=np.ndarray>

尝试用以下代码绘制时，输出为近乎方形的灰色图像，未达到预期效果：

import matplotlib.pyplot as plt

# Squeeze the dimensions of the satellite image
image = df['data'].squeeze()

# Plot the satellite image
plt.imshow(image, origin='lower')
plt.show()

问题分析 #

1. Dask延迟计算未触发：data是Dask数组，squeeze后仅保留元数据，未加载实际像素值，导致显示异常。
2. 图像比例失真：原始图像为4176×5548的矩形，imshow默认等比例显示，受窗口限制被压缩成方形。
3. 色阶未适配卫星数据：HRV数据的动态范围未被正确映射，默认色阶导致整体偏灰。
4. 地理坐标未利用：直接用imshow忽略了静止卫星的投影坐标，无法准确反映图像的地理范围。

解决步骤 #

1. 加载实际数据（处理Dask延迟计算） #

Dask数组需要触发计算才能获取实际像素值，可完整加载或先取子集测试：

# 完整加载数据（注意：大尺寸数据需确保内存足够）
image_data = df['data'].squeeze().compute()

# 子集测试（避免内存占用过高）
# image_data = df['data'].isel(y_geostationary=slice(0,1000), x_geostationary=slice(0,1000)).squeeze().compute()

2. 校正图像显示比例 #

通过aspect='auto'保持原始宽高比，或结合地理坐标设置显示范围：

import matplotlib.pyplot as plt

plt.figure(figsize=(12, 8))  # 调整画布大小适配图像比例

# 方法1：仅保持原始比例
plt.imshow(image_data, origin='lower', aspect='auto')

# 方法2：结合投影坐标设置extent（更精准）
x_min = df['x_geostationary'].min().compute()
x_max = df['x_geostationary'].max().compute()
y_min = df['y_geostationary'].min().compute()
y_max = df['y_geostationary'].max().compute()
plt.imshow(image_data, origin='lower', extent=[x_min, x_max, y_min, y_max], aspect='auto')

3. 适配HRV数据的色阶 #

HRV数据需调整色阶避免极端值影响，可基于分位数设置vmin和vmax：

# 查看数据统计值，确定有效范围
print(f"数据范围：{image_data.min()} ~ {image_data.max()}，均值：{image_data.mean()}")

# 基于2%和98%分位数设置色阶，过滤极端值
vmin = image_data.quantile(0.02)
vmax = image_data.quantile(0.98)

plt.imshow(image_data, origin='lower', aspect='auto', vmin=vmin, vmax=vmax, cmap='Greys_r')
plt.colorbar(label='HRV反射率')  # 添加色标解释数值含义

4. 使用xarray内置plot方法（更简便） #

xarray自带的plot.imshow可自动处理维度和坐标，无需手动操作：

df['data'].squeeze().plot.imshow(
    origin='lower',
    aspect='auto',
    vmin=vmin,
    vmax=vmax,
    cmap='Greys_r',
    figsize=(12,8)
)
plt.title('HRV卫星图像 - 2021-01-01')
plt.show()

5. 大尺寸数据优化方案 #

若数据过大无法完整加载，可直接用Dask数组渲染，或结合地理投影工具：

• 直接调用Dask数组的plot方法：df['data'].squeeze().plot.imshow(...)（无需提前compute）
• 使用cartopy添加地理底图，适配静止卫星投影

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