spatial_spectral_smooth.py

import numpy as np
from astropy import units as u
from astropy.convolution import Gaussian1DKernel, Gaussian2DKernel
import os
import glob
from spectral_cube import SpectralCube


for fn in [
    "G001_13CO21_SEDIGISM_DR1c.fits",
    "G002_13CO21_SEDIGISM_DR1c.fits",
    "G003_13CO21_SEDIGISM_DR1c.fits",
    "G004_13CO21_SEDIGISM_DR1c.fits",
    "G005_13CO21_SEDIGISM_DR1c.fits",
    "G359_13CO21_SEDIGISM_DR1c.fits",
    "G358_13CO21_SEDIGISM_DR1c.fits",
    "G357_13CO21_SEDIGISM_DR1c.fits",
    "G356_13CO21_SEDIGISM_DR1c.fits",
    "G355_13CO21_SEDIGISM_DR1c.fits",
    "G000_13CO21_SEDIGISM_DR1c.fits",
]:
    for smfactor in (2,4,8):
        outfn = fn.replace(".fits", "_spatialdownsample_{0}x.fits".format(smfactor))
        if not os.path.exists(outfn):

            print(fn, outfn)

            cube = SpectralCube.read(fn)

            fwhm_factor = np.sqrt(8*np.log(2))
            kernel = Gaussian2DKernel(smfactor)

            sm_cube = cube.spatial_smooth(kernel)
            ds_cube = sm_cube[:,::smfactor,::smfactor]

            ds_cube.write(outfn)

    dsfn = fn.replace(".fits", "_spatialdownsample_4x.fits")

    for target_resolution in (1,2,5):
        outfn = dsfn.replace(".fits", "_downsampled_{0}kms.fits".format(target_resolution))
        if not os.path.exists(outfn):

            print(fn, dsfn, outfn)

            cube = SpectralCube.read(dsfn).with_spectral_unit(u.km/u.s)

            fwhm_factor = np.sqrt(8*np.log(2))
            hanning_factor = 1129/977.
            current_resolution = np.mean(np.diff(cube.spectral_axis)) * hanning_factor
            target_resolution = u.Quantity(target_resolution, u.km/u.s)
            pixel_scale = current_resolution
            gaussian_width = ((target_resolution**2 - current_resolution**2)**0.5 /
                              pixel_scale / fwhm_factor)
            kernel = Gaussian1DKernel(gaussian_width)

            new_xaxis = np.arange(-200, 200, target_resolution.value/2) * u.km/u.s
            print("Beginning spectral_smooth")
            sm_cube = cube.spectral_smooth(kernel, num_cores=8)
            print("Beginning spectral_interpolate")
            ds_cube = sm_cube.spectral_interpolate(new_xaxis)

            ds_cube.write(outfn)