code: def extract_mapcube_region(infile, skydir, outfile, maphdu=0):
      
  
      h = fits.open(os.path.expandvars(infile))
  
      npix = 200
      shape = list(h[maphdu].data.shape)
      shape[1] = 200
      shape[2] = 200
  
      wcs = WCS(h[maphdu].header)
      skywcs = WCS(h[maphdu].header, naxis=[1, 2])
      coordsys = get_coordsys(skywcs)
  
      region_wcs = wcs.deepcopy()
  
      if coordsys == 'CEL':
          region_wcs.wcs.crval[0] = skydir.ra.deg
          region_wcs.wcs.crval[1] = skydir.dec.deg
      elif coordsys == 'GAL':
          region_wcs.wcs.crval[0] = skydir.galactic.l.deg
          region_wcs.wcs.crval[1] = skydir.galactic.b.deg
      else:
          raise Exception('Unrecognized coordinate system.')
  
      region_wcs.wcs.crpix[0] = npix // 2 + 0.5
      region_wcs.wcs.crpix[1] = npix // 2 + 0.5
  
      from reproject import reproject_interp
      data, footprint = reproject_interp(h, region_wcs.to_header(),
                                         hdu_in=maphdu,
                                         shape_out=shape)
  
      hdu_image = fits.PrimaryHDU(data, header=region_wcs.to_header())
      hdulist = fits.HDUList([hdu_image, h['ENERGIES']])
      hdulist.writeto(outfile, clobber=True)
docs: Extract a region out of an all-sky mapcube file.
  
      Parameters
      ----------
  
      infile : str
          Path to mapcube file.
  
      skydir : `~astropy.coordinates.SkyCoord`