pyspeckit / pyspeckit

Python Spectroscopic Toolkit
MIT License
104 stars 83 forks source link

Provide minimal example of cube fitting #157

Open astrofrog opened 8 years ago

astrofrog commented 8 years ago

I've had a couple of users comment that it's not easy to fit all pixels in spectral cubes. I pointed them to some of the examples, but these are pretty complex, e.g:

http://pyspeckit.readthedocs.org/en/latest/example_nh3_cube.html

It would be nice to provide an example that shows the minimal amount of code required to fit a single cube, since in these long examples it isn't easy to see what's optional and what's required.

keflavich commented 8 years ago

This one is simpler: http://pyspeckit.readthedocs.org/en/latest/example_cube_fromscratch.html

Would it be better to separate out the "make a cube" and "fit a cube" sections?

astrofrog commented 8 years ago

@keflavich - yes, I think that would help!

keflavich commented 8 years ago

OK, it's on the todo list. Have you by any chance heard from advanced users who might be interested in contributing some examples?

astrofrog commented 8 years ago

@keflavich - no, but I'll be on the lookout

vlas-sokolov commented 8 years ago

Would this do? It's small but should still be useful for new users:

import pyspeckit
from pyspeckit.cubes.tests.test_cubetools import make_test_cube
import matplotlib.pylab as plt
import numpy as np

# generate a test spectral cube (10x10, with a 100 spectral channels)
make_test_cube((100,10,10), outfile='test.fits')
spc = pyspeckit.Cube('test.fits')

# do a crude noise estimate on the 30 edge channels
rmsmap = np.vstack([spc.cube[:15], spc.cube[85:]]).std(axis=0)

# get a cube of moments
spc.momenteach(vheight=False)

# fit each pixel taking its moment as an initial guess
spc.fiteach(fittype = 'gaussian',
            guesses = spc.momentcube,
            errmap = rmsmap,
            signal_cut = 3, # ignore pixels with SNR<3
            blank_value = np.nan,
            start_from_point=(5,5))
spc.mapplot()
# show the fitted amplitude
spc.show_fit_param(0, cmap='viridis')
plt.show()
keflavich commented 8 years ago

...yes! that's good!

...seriously, there's a show_fit_param method? I have completely forgotten about this. I had also completely forgotten there is a make_test_cube method. Thanks Vlas! Make this a PR!