import regions
import glob
from astropy import units as u
from astropy.coordinates import SkyCoord
from spectral_cube import SpectralCube

center = SkyCoord('17:47:24.4487906862 -28:21:43.5028185481', unit=(u.h, u.deg))
reg = regions.RectangleSkyRegion(center=center, width=35*u.arcsec, height=35*u.arcsec)

for fn in glob.glob("full_SgrB2_TETC7m*lines.fits"):
    print(fn)
    cube = SpectralCube.read(fn)
    mcube = cube.subcube_from_regions([reg])
    mcube.write(fn.replace(".fits","_busch_q.fits"),
                overwrite=True)