I'd like to plot an array with the following dimensions: [3x3x180x360]
it's a 180x360 world map with 3 different opacity layers and 3 different pressure levels. Thus, I would like to plot the map with data of 1 opacity layer at one pressure level.
Using the ReadAsArray() function of gdal, results in loosing one dimension: [64800, 3, 3]
Is there another way or an argument to keep the 4 dimensions?
import numpy as np from osgeo import gdal 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).ReadAsArray() # takes ~5min to read ... #filtering bad data: cld_amount_liq_md[cld_amount_liq_md > 3.40E38] = np.nan # Plotting plt.imshow(cld_amount_liq_md ,cmap ="jet")```