code: def extract_images_from_tscube(infile, outfile):
      
      inhdulist = fits.open(infile)
      wcs = pywcs.WCS(inhdulist[0].header)
      map_shape = inhdulist[0].data.shape
  
      t_eng = Table.read(infile, "EBOUNDS")
      t_scan = Table.read(infile, "SCANDATA")
      t_fit = Table.read(infile, "FITDATA")
  
      n_ebin = len(t_eng)
      energies = np.ndarray((n_ebin + 1))
      energies[0:-1] = t_eng["E_MIN"]
      energies[-1] = t_eng["E_MAX"][-1]
  
      cube_shape = (n_ebin, map_shape[1], map_shape[0])
  
      wcs_cube = wcs_utils.wcs_add_energy_axis(wcs, energies)
  
      outhdulist = [inhdulist[0], inhdulist["EBOUNDS"]]
  
      FIT_COLNAMES = ['FIT_TS', 'FIT_STATUS', 'FIT_NORM',
                      'FIT_NORM_ERR', 'FIT_NORM_ERRP', 'FIT_NORM_ERRN']
      SCAN_COLNAMES = ['TS', 'BIN_STATUS', 'NORM', 'NORM_UL',
                       'NORM_ERR', 'NORM_ERRP', 'NORM_ERRN', 'LOGLIKE']
  
      for c in FIT_COLNAMES:
          data = t_fit[c].data.reshape(map_shape)
          hdu = fits.ImageHDU(data, wcs.to_header(), name=c)
          outhdulist.append(hdu)
  
      for c in SCAN_COLNAMES:
          data = t_scan[c].data.swapaxes(0, 1).reshape(cube_shape)
          hdu = fits.ImageHDU(data, wcs_cube.to_header(), name=c)
          outhdulist.append(hdu)
  
      hdulist = fits.HDUList(outhdulist)
      hdulist.writeto(outfile, clobber=True)
      return hdulist
docs: Extract data from table HDUs in TSCube file and convert them to FITS images