CASA user's guide to common operations in spectral-cube

[e-koch]

Split out tutorial initially in #11. Depends on pending PR radio-astro-tools/spectral-cube#665.

Temporary binder link:

Link to main repo is in the README.

[e-koch]

The install on binder is failing. Will need to revise the build files

[e-koch]

Build still not working. Seem like:

• Building spectral-cube and radio-beam leads pip to also try to build casatasks and casatools, that then fails on finding setupstools_scm.
• With only casatools/casatasks, the environment builds but neither can be imported when the notebook server starts

[e-koch]

Yay! CASA is finally working on binder!

Here's what ended up being required for the install:

• apt.txt to install libgfortran3
• runtime.txt to specify python 3.6 should be installed (current default right now is 3.7)
• requirements.txt List the casa pip mirror with --extra-index-url. Otherwise using --index-url with the casa mirror causes the install of spectral-cube and radio-beam fail because that mirror lacks a few requirements

Other things to note:

• the build of pyregion seems to fail. This looks like it might be a numpy version clash but I didn't copy out the build log for the current working binder image before it launched the notebook

[e-koch]

Ready once radio-astro-tools/spectral-cube#680 is merged.

[keflavich]

ran locally, and... something's amiss?

Amplitude: spectral-cube: 12.871325161920732 K km / s, CASA: 12.8383888569787 K km / s
Major stddev: spectral-cube: 4.927549702771874, CASA: 4.451009139383827
Minor stddev: spectral-cube: 3.6335649764546365, CASA: 3.4764324459573515
Position angle: spectral-cube: -26.70443925746928 deg, CASA: 132.2080778767417 deg

The Gaussian shapes are different, but they look like they evaluated to be the same

[e-koch]

Yeah, they're within ~1/5 of a pixel (1 pixel = 1" here). I wasn't able to get these to agree any better, but I agree that they seem to be evaluated almost the same.

I can add a note about this in the notebook.

[keflavich]

[keflavich]

I think these are too different. I'm concerned now.

[keflavich]

The difference is with the underlying image being fit. CASA seems to only be able to cutout odd-sized chunks and is dealing with masking differently.

[keflavich]

the odd-size chunk issue can be solved with box instead of centerbox:
region=f'box [ [{xcasa-size}pix,{ycasa-size}pix] , [{xcasa+size-1}pix, {ycasa+size-1}pix]]', # Slice out location of source

The remaining difference might be that astropy masks pixels (ignores them), casa doesn't

[keflavich]

nope, my idea is wrong. You've manually set bad pixels to zero. So what is CASA doing?

[keflavich]

the difference is in the fitter, verified by astropy-fitting the CASA data:

ia.open(outfilename)
casa_mom0 = ia.getchunk().squeeze().T
ia.close()
p_gauss2D = fit_p(p_init_gauss2D, xx, yy, casa_mom0*u.K*u.km/u.s)

it gives the same result as astropy-fitting the spectral-cube data

[keflavich]

This MWE works as expected:

import numpy as np
import shutil
import os
from casatools import image; ia = image()
from casatasks import imfit
from astropy.modeling.models import Gaussian2D
from astropy.modeling.fitting import LevMarLSQFitter
from astropy import units as u
FWHM2SIGMA = (8*np.log(2))**0.5
mod = Gaussian2D(amplitude=1.5, x_mean=15, y_mean=17, x_stddev=3.2, y_stddev=1.5, theta=45*u.deg)
yy,xx = np.indices([64,64])
modim = mod(xx,yy)
if os.path.exists('test_gauss.image'):
    shutil.rmtree('test_gauss.image')
ia.fromarray(outfile='test_gauss.image', pixels=modim.T)
spatfit_casa = imfit('test_gauss.image')
g_casa_spatfit = Gaussian2D(amplitude=spatfit_casa['results']['component0']['peak']['value'],
                                   x_mean=spatfit_casa['results']['component0']['pixelcoords'][0],
                                   y_mean=spatfit_casa['results']['component0']['pixelcoords'][1],
                                   y_stddev=spatfit_casa['results']['component0']['shape']['majoraxis']['value'] / FWHM2SIGMA / 60,
                                   x_stddev=spatfit_casa['results']['component0']['shape']['minoraxis']['value'] / FWHM2SIGMA / 60,
                                   theta=spatfit_casa['results']['component0']['shape']['positionangle']['value'] * u.deg
                                   )
print(g_casa_spatfit, mod)

[keflavich]

MWE with noise:

import numpy as np
import shutil
import os
from casatools import image; ia = image()
from casatasks import imfit
from astropy.modeling.models import Gaussian2D
from astropy.modeling.fitting import LevMarLSQFitter
from astropy import units as u
FWHM2SIGMA = (8*np.log(2))**0.5
mod = Gaussian2D(amplitude=1.5, x_mean=15, y_mean=17, x_stddev=3.2, y_stddev=1.5, theta=45*u.deg)
yy,xx = np.indices([64,64])

np.random.seed(12345)
noise = np.random.randn(64,64) / 2
modim = mod(xx,yy) + noise
if os.path.exists('test_gauss.image'):
    shutil.rmtree('test_gauss.image')
ia.fromarray(outfile='test_gauss.image', pixels=modim.T)
ia.close()
spatfit_casa = imfit('test_gauss.image')
g_casa_spatfit = Gaussian2D(amplitude=spatfit_casa['results']['component0']['peak']['value'],
                                   x_mean=spatfit_casa['results']['component0']['pixelcoords'][0],
                                   y_mean=spatfit_casa['results']['component0']['pixelcoords'][1],
                                   y_stddev=spatfit_casa['results']['component0']['shape']['majoraxis']['value'] / FWHM2SIGMA / 60,
                                   x_stddev=spatfit_casa['results']['component0']['shape']['minoraxis']['value'] / FWHM2SIGMA / 60,
                                   theta=spatfit_casa['results']['component0']['shape']['positionangle']['value'] * u.deg
                                   )
fitter = LevMarLSQFitter()
pguess_astropy = Gaussian2D(amplitude=2, x_mean=16, y_mean=16, x_stddev=2, y_stddev=2, theta=0*u.deg)
pfit_astropy = fitter(pguess_astropy, xx, yy, modim)
print(f"CASA: {g_casa_spatfit}\n\n Astropy: {pfit_astropy}\n\n input:{mod}")

[e-koch]

@keflavich -- thanks for digging into the fitting problem.

Also not due to imfit doing something with the beam size. Deleting the beam in the CASA moment 0 image makes no difference.

    ia.open(outfilename)
    ia.setrestoringbeam(remove=True)
    print(ia.restoringbeam())
    ia.close()

[e-koch]

Also no difference for adding the same initial guesses to the CASA fit that we're using for the astropy fit:

# Create a new moment0 image
outfilename = f"{filename.split('.')[0]}.moment0.image"

# Cutout a 13 pixel box around the central source
size = 13

# Identify the spectral channels. We will only integrate over the range where
# there is bright signal
low_chan = cube.closest_spectral_channel(-230 * u.km / u.s)
high_chan = cube.closest_spectral_channel(-190 * u.km / u.s)

os.system(f"rm -r {outfilename}")

print(x, y)

if not os.path.exists(outfilename):
    print('Making new moment image')
    # Slice out the same size we will use with spectral-cube
    # See the region file formats here: https://casaguides.nrao.edu/index.php/CASA_Region_Format
    immoments(filename, moments=0, axis='spec',
              outfile=outfilename,
              mask=f'"{filename}">0.12',  # Mask out noise
              chans=f'{low_chan}~{high_chan}',
#               region=f'box [ [{x-size}pix,{y-size}pix] , [{x+size-1}pix, {y+size-1}pix]]')
              region=f'centerbox [ [{x}pix,{y}pix] , [{2*size}pix, {2*size}pix]]', # Slice out location of source
              )

    ia.open(outfilename)
    ia.setrestoringbeam(remove=True)
    print(ia.restoringbeam())
    ia.close()

# Make a text file with the guesses
if os.path.exists('casafitguesses.txt'):
    os.system('rm casafitguesses.txt')
os.system("touch casafitguesses.txt")
os.system("echo -n '12, 13, 13, 8arcsec, 8arcsec, 0deg' >> casafitguesses.txt")

with open('casafitguesses.txt') as f:
    print(f.readlines())


# Fit 1 2D Gaussian to the cropped moment0 image made above.
spatfit_casa = imfit(imagename=outfilename, estimates='casafitguesses.txt')

# Note that FWHM needs to be converted to "sigma" (FWHM = sqrt(8 ln 2) sigma).
FWHM2SIGMA = np.sqrt(8 * np.log(2))

# Make a grid of spatial grids in terms of the pixel shape. We're going to leave out the units for now.
yy_pix, xx_pix = np.mgrid[:2 * size, :2 * size]

t0 = time()

# There is some discrepancy in units in the format that CASA returns.
# For this example, we will only use the peak amplitude unit K (though CASA returns this as Jy/pixel)
g_casa_spatfit = models.Gaussian2D(amplitude=spatfit_casa['results']['component0']['peak']['value'] * u.K * u.km / u.s,
                                   x_mean=spatfit_casa['results']['component0']['pixelcoords'][0],
                                   y_mean=spatfit_casa['results']['component0']['pixelcoords'][1],
                                   y_stddev=spatfit_casa['results']['component0']['shape']['majoraxis']['value'] / FWHM2SIGMA,
                                   x_stddev=spatfit_casa['results']['component0']['shape']['minoraxis']['value'] / FWHM2SIGMA,
                                   theta=spatfit_casa['results']['component0']['shape']['positionangle']['value'] * u.deg)

t1 = time()

print(f"CASA imfit time: {t1 - t0} s")


print("CASA 2D Gaussian fit:")
print(g_casa_spatfit)

[keflavich]

Can it be something to do with the masks? I tried masking out the data in the astropy fit (your version doesn't mask out the data - suggest modifying so that it does; set weights=~np.isnan(data) before setting nans to zero), but that didn't make them match either. I didn't get around to testing whether the mask in the casa imfit is causing this. Maybe CASA masks out any spatial pixel that has any masked out volume pixel?

[e-koch]

Closing the loop on the fit discrepancies and @keflavich's MWEs:

• the issue was in the astropy fitting when fitting to the WCS coordinate grid and handling the cos dec correction. It's a lot easier to just fit to the pixel grid, then convert pixel to angular scales
• There was no noticeable difference in changing the masking/weights. So it's unclear if there is a difference in the mask handling for CASA vs. astropy.modeling, but it doesn't lead to any significant differences

I've also cleaned up the box selections in CASA to use box instead of centerbox. There was a +1 offset in the spectral fit selection, too, that's all fixed now.

@keflavich -- Mind having one last look through the spatial fitting? The key addition is handling the pix -> ang scale conversion. There's a bit of confusion with the y_stddev/x_stddev naming matching up with the y, x pixel scales that I've tried clarifying

[e-koch]

@keflavich Good to squash and merge (the first few commits included the cell outputs in the nb)