How can I overlay an image (calculated from arrays) on the map (Mercator proj) by coordinates?

I have a HDF file with 3 arrays for the image, c 2 for the coordinates (lat lon). Arrays of size 2784 2784.

I'm trying to draw an RGB image:

import cartopy  import cartopy.crs as ccrs  import skimage.transform as st  import matplotlib.pyplot as plt  import matplotlib.cm as mcm  from matplotlib.colors import ListedColormap  import h5py  import numpy as np    def test(hdfFile):  # latitude, longitude, image      satChan = ['B2_0.72', 'BT4_3.75', 'BT9_10.7']      lat = hdfFile['BT_Lat'][:, :]      lon = hdfFile['BT_Lon'][:, :]      satChanLen = len(satChan)      dafa = []      mask = np.zeros((2784, 2784))      for Nch in range(satChanLen):          ds = satChan[Nch]          data_set = hdfFile[ds]          if ds == 'B2_0.72':              data = st.resize(data_set[:, :], (2784, 2784))              mask = (data < 0.5) * (data > 0)  # & data>0              # data[mask] = 1              dafa.append(data)          else:              data = data_set[:, :]              dafa.append(data)      daf = np.array(dafa)      porT = 0.03  # переход к ассимптоте      y1 = porT * 255      y2 = (1 - porT) * 255      k = 255      lo = np.log(y1 / (1 - y1 / k)) / np.log(y2 / (1 - y2 / k))      lu = np.log(y1 / (1 - y1 / k))      v = daf.shape      x, y = v[1], v[2]      arr = np.zeros((x, y, 3))      for Nch in range(satChanLen):          masked = np.ma.masked_equal(daf[Nch, :, :], 0)          ssi = masked.min()          ssa = masked.max()          step1 = np.where(daf[Nch, :, :] > 0, 255 * (daf[Nch, :, :] - ssi) / (ssa - ssi), 0)            histShag = 10          porog = 0.1          b = np.arange(0, 300, histShag)            histStep = step1.flat[step1.flat > 0]          th = np.histogram(histStep, b)          por = max(th[0]) * porog          mm = []          for i in range(int(300 / histShag) - 1):              if th[0][i] > por:                  mm.append(th[1][i])          mmmin = min(mm)          mmmax = max(mm)          ak = (mmmin - mmmax * lo) / (1 - lo)          sk = lo * (mmmax - mmmin) / (lu * (1 - lo))          exp = np.where(step1 > 0, np.exp((step1 - ak) / sk), 0)          # step2 = numpy.where(step1 > 0, exp / (1 + exp / 254), 0)          # step2 = numpy.where(step2 == 0, 255, step2)          if Nch < 2:              step2 = np.where(step1 > 0, exp / (1 + exp / 254), 0)          else:              step2 = np.where(step1 > 0, exp / (1 + exp / 255), 0)          print(step2)          step2[mask] = 128          arr[:, :, Nch] = step2      arr = arr.astype("uint8")        extent = (-60, 80, 30, 90)  # Left Longitude, Right longitude, Down latitude, Up latitude      proj = ccrs.Mercator(central_longitude=40)      ax = plt.axes(projection=proj)      # ax.add_feature(cartopy.feature.LAND)      ax.add_feature(cartopy.feature.OCEAN)      ax.add_feature(cartopy.feature.COASTLINE)      ax.add_feature(cartopy.feature.BORDERS, linestyle='-', alpha=.5)      ax.add_feature(cartopy.feature.LAKES, alpha=0.95)      ax.add_feature(cartopy.feature.RIVERS)      ax.set_extent(extent, crs=ccrs.PlateCarree())      ax.gridlines(draw_labels=True)      # ax.pcolormesh(lon, lat, c=rgb, transform=ccrs.PlateCarree())      ax.scatter(lat, lon, c=arr[0], transform=ccrs.PlateCarree())      #ax.contourf(lat, lon, arr, 1, transform=ccrs.PlateCarree())      plt.show()    if __name__ == '__main__':      path = r"202104230100_3720.h5"      with h5py.File(path, 'r') as f:          test(f)          # lat = f["BT_Lat"][:]          # lon = f["BT_Lon"][:]          # img = f["BT7_8.70"][:]          # print(lat.shape)          # print(lon.shape)          # print(img.shape)  

I get exceptions, I tried all the options (contourf, pcolormesh, scatter).

How do I overlay an RGB image on a grid?

Test code, for drawing 1 channel. It works if you draw on 1 channel. I tried to distill all 3 channels in list((r,g,b)) element by element, but it also does not work.

TestCode:

def test2(hdf):      satChan = ['B2_0.72', 'BT4_3.75', 'BT9_10.7']      lat = hdf['BT_Lat'][:, :]      lon = hdf['BT_Lon'][:, :]      satChanLen = len(satChan)      x, y = lat.shape      dafa = np.zeros([x, y, 3], dtype=np.uint8)      for Nch in range(satChanLen):          ds = satChan[Nch]          data_set = hdf[ds]          if ds == 'B2_0.72':              dafa[:, :, Nch] = st.resize(data_set[:, :], (2784, 2784))          else:              dafa[:, :, Nch] = data_set[:, :]      # extent = (-60, 80, 30, 90)  # Left Longitude, Right longitude, Down latitude, Up latitude      extent = (lon[lon > -999].min(), lon.max(), lat[lat > -999].min(), lat.max())      proj = ccrs.Mercator(central_longitude=40)      ax = plt.axes(projection=proj)      # ax.add_feature(cartopy.feature.LAND)      ax.add_feature(cartopy.feature.OCEAN)      ax.add_feature(cartopy.feature.COASTLINE)      ax.add_feature(cartopy.feature.BORDERS, linestyle='-', alpha=.5)      ax.add_feature(cartopy.feature.LAKES, alpha=0.95)      ax.add_feature(cartopy.feature.RIVERS)      ax.set_extent(extent, crs=ccrs.PlateCarree())      ax.gridlines(draw_labels=True)      ax.scatter(lon, lat, c=dafa[:, :, 2], transform=ccrs.PlateCarree())      plt.show()  

File: Source file(hdf) ~ 3Гб

https://stackoverflow.com/questions/67223918/how-can-i-overlay-an-image-calculated-from-arrays-on-the-map-mercator-proj-b April 23, 2021 at 12:14PM