Python GDAL读取4D数组的索引选取及plt.imshow绘图实现

• 待可视化目标数组逻辑维度为[3×3×180×360]，对应分辨率180×360的全球云量地图，包含3个不透明度层、3个气压层两个分类维度。需求为选取单个不透明度层、单个气压层对应的二维地图数据，通过plt.imshow()完成可视化，Panoply软件下的该数据可视化效果可参考附图。

import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt

granule = "CER_CldTypHist_GEO-MODIS_Edition4A_407408.202109.hdf"
hdf_file = gdal.Open(workdir_data + "/" + granule)
subDatasets = hdf_file.GetSubDatasets()
cld_amount_liq_md = gdal.Open(subDatasets[68][0]).ReadAsArray() # 读取耗时约5分钟

# ReadAsArray()读取后数组形状为[64800, 3, 3]，因此执行重整形：
cld_amount_liq_md  = cld_amount_liq_md.reshape(180,360, 3, 3)

# 过滤无效值：
cld_amount_liq_md[cld_amount_liq_md > 3.40E38] = np.nan

# 绘图环节：如何正确索引数组完成绘图？
plt.imshow(cld_amount_liq_md[?,?,?] ,cmap ="jet")

reshape得到的形状为(180,360,3,3)的四维数组，应如何填写索引选取指定气压层、指定不透明度层的二维切片，传入plt.imshow()完成全球云量地图绘制。

• 可通过以下代码遍历获取HDF文件的子数据集维度等元信息，确认维度对应关系：

from osgeo import gdal

granule = "CER_CldTypHist_GEO-MODIS_Edition4A_407408.202109.hdf"
hdf_file = gdal.Open(workdir_data + "/" + granule)
subDatasets = hdf_file.GetSubDatasets()

j=0
for i in subDatasets:
    print(j, i[1]) # i[0]为子数据集路径，i[1]为子数据集元信息
    j+=1

• 目标子数据集全变量名为：Monthly_Day_Averages/Cloud_Properties_for_9_Cloud_Types_Monthly_Day/cld_amount_liq_md

当前reshape得到的四维数组维度顺序依次为纬度（180格点）、经度（360格点）、不透明度层（3级）、气压层（3级），索引规则如下：

1. 前两个经纬度维度需要全选，用:表示保留所有格点
2. 后两个维度填入目标层对应的索引值，索引从0开始计数，对应第1到第3级分类
3. 如果绘图后发现层对应关系和Panoply不一致，调换后两个维度的索引顺序即可（即气压层在前、不透明度层在后），可结合元信息遍历输出的维度说明确认顺序。

完整绘图示例代码：

# 按需修改两个层的索引，取值范围0、1、2，对应3个分级
opacity_level_idx = 0  # 不透明度层索引
pressure_level_idx = 0 # 气压层索引

# 提取二维切片
cld_2d = cld_amount_liq_md[:, :, opacity_level_idx, pressure_level_idx]

# 绘图：origin='lower'将坐标原点设为左下角，匹配地理坐标从南到北的顺序；extent指定经纬度范围
plt.figure(figsize=(12,6))
plt.imshow(cld_2d, cmap="jet", origin="lower", extent=[-180, 180, -90, 90])
plt.colorbar(label='液态水云量')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

若切片提取后出现东西半球错位、纬度反向的问题，可通过np.flip()、np.roll()对二维数组做方向调整，匹配Panoply的显示效果即可。

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