1

I have a fits datacube with galactic longitude, latitude and velocity in the 3 axis. To extract the spectrum from the datacube at a particular pixel value of longitude and latitude, I use the function

cube[:, 1935, 1407].quicklook()
plt.show()

and the image is extracted with the function

cube.to_pvextractor()
plt.show()

A sample spectrum spectrum and a zoomed image image is attached here.

The bright spots are the detections. How do I use several pixels and average the spectra to get a mean spectrum so that I reduce the noise and analyze the peak? I have been trying to code this but I don't know how to proceed as I am new to python. Can anybody please give a hint?

Lidia
  • 27
  • 5
  • How does upper image relate to lower-left? Why lower-right is blank? Also: the background data in the lower-left has significant noise from the instruments (clear diagonal "strips") that has to be accounted for. The most important thing: astronomical data should be processed the way the astronomers do. I'm afraid SO is not the best place to ask this question. – zkoza May 25 '21 at 11:32
  • Thank you for pointing out the problem in my question. Just an additional query. Is there any way to take the average of several pixels and form a new spectrum? – Lidia May 25 '21 at 13:06
  • I'm not an astronomer. If I were to tackle this all by myself, I'd do it as in 1D case: fit a 2D Gaussian to a noisy 2D signal. Actually, there are at least 4 distinct signal sources (galaxies?) in the data. You can pinpoint them one after another, eventually you'll have a superposition of 4 Gaussians. You need python's scipy. Then manually estimate the result (coordinates of the signal center, its "width" and "amplitude". Run sci-py. Take your PhD :-) BTW, are you a PhD student or is it only Master's in astronomy? – zkoza May 25 '21 at 14:30
  • I am doing a short-term project in astronomy after completing my M.Sc in Physics. I am new to python and you have been a great help to me. Thank you very much. – Lidia May 25 '21 at 15:55
  • 1
    Then remember to estimate the initial values of the fitted parameters, otherwise sci-py may get lost and return a completely absurd solution. There's nothing wrong in this. When sci-py returns the solution, you must compare it against the original, experimental data to make sure sci-py has converged to the required solution. – zkoza May 25 '21 at 17:33
  • Suppose if I create a `subcube = cube[:, 1450:1620, 430:500]` such that it includes the bright spot and the write `spectrum = subcube.mean(axis=(1,2))` will it give me an averaged spectrum? – Lidia May 26 '21 at 06:58
  • 1
    Yes and no. Yes, means, this would be a very crude estimate, most likely a highly underestimated value. No, because fitting data to a 2D Gaussian is rather simple and robust. sci-py has all you need. Unfortunately, I don't use python nor sci-py, but I've seen it many times in action, it is a full replacement of octave/matlab, so nothing more is needed. – zkoza May 26 '21 at 14:29

1 Answers1

0

You can use spectrapepper for this:

import spectrapepper as spep
import matplotlib.pyplot as plt

# load sample data from library
x, y = spep.load_spectras()

# calculate the average spectra of the set
avg = spep.avg(y)

# plot the result compared to the data set
for i in y:
    plt.plot(x, i, lw=0.5)
plt.plot(x, avg, c='red', lw=1)
plt.show()

The library has also other tools for data set analysis. Check example 6 included in the docs for more options of similar nature.

DharmanBot
  • 1,066
  • 2
  • 6
  • 10
Enric Grau-Luque
  • 114
  • 11