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[68][0]).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")```

1 Answer 1


I found the solution here: https://gis.stackexchange.com/a/315814

you have to reshape the array with band.ReadAsArray().reshape(3,3,180,360)

Not the answer you're looking for? Browse other questions tagged or ask your own question.