Q3D / q3dfit

Python software for fitting integral field spectroscopic data
GNU General Public License v3.0
10 stars 3 forks source link

q3dfit NIRSpec spectrum extraction #11

Closed WujiWang-astro closed 1 year ago

WujiWang-astro commented 1 year ago

Hi, I am using q3dfit on NIRSpec data cube and simply following the notebook. On the first try, I found that there are several parts where more description could be provided:

drupke commented 1 year ago

Hi Wuji,

Thanks so much for these questions.

what can we specify in the q3di.argsreadcube? This is a dictionary of arguments that can be passed to the readcube.Cube class, which reads in the data cube when it is instantiated. The load_cube() method of the q3din class passes these on when it is invoked. So any a related question: how does the flux conversion work? The flux conversion is a work in progress, and I agree not fully transparent. By default, Cube expects input fluxes in MJy/sr (the default output of the JWST pipeline). The code (except when using the questfit continuum fitter — more on that to come in future releases) then converts the input to, and expresses any output in, units of erg/s/cm^2/micron/sr. (Equivalently, both the input /sr and output /sr can be ignored.) For example, when I specify both the 'fluxunit_in' and 'fluxunit_out' as 'erg/s/cm2/Angstrom', do I expect there is no conversion going on? But in this case the output flux level is then several orders of magnitude different than what I have in the cube. This appears to be a bug. I am getting ready to issue a patch release in the next week or so with a few changes, so I will make this fix and include it. Thanks for finding this. then this leads to the third question: the spectra extracted by 'cube.specextract' is at which location? The description of the col and row is not clear. Also, the unit of the radius given is not clear either in the documentation. Does it extract a spectrum at the spaxel location of (col, row) if we give radius=0? I agree that this is a bit unclear. The first two arguments are the (unity-offset) spaxel location at which the spectrum is extracted using a circular aperture. The radius of the aperture, again in spaxels, is given by the third argument. The extraction is done by default using aperture photometry with photutils.aperture.aperture_photometry(). Here you have identified a second bug, however. The radius=0 aperture then defaults to radius=0.01, which using the default “exact” calculation from aperture_photometry() results in a literal aperture of this size, not simply a single-spaxel extraction. I’ll fix this bug in the upcoming patch release, as well.

Dave

WujiWang-astro commented 1 year ago

Hi Dave, Thank you very much for the reply. q3dfit is indeed a nice and convenient tool for fitting. Looking forward to the new release. Another issue which I am not sure is universal or just a local problem on my laptop: I followed the example of nirspec notebook and quickly ran through the fitting. In general, everything works fine for me despite some untransparent details. However, when I call the plot_line function, it returns me the error that variable 'fcn' is not str (in line 628 of q3dout.py). This is wired because it is specified to be one in line 618. Right now I specify it by myself again, fcn='plotline', and solve it. But it should be worth posting it in case something should be fixed. Wuji

WujiWang-astro commented 1 year ago

Update of the cube.specextract: when I give a radius of 3 pix, and specify the fluxunit_in and fluxunit_out as the same, then the result is on the consistent level as I expect. So I think the previous problem was just an issue of the 'radius=0' (or 0.01) thing. In general this works fine, just need a little bit more description of how this works.

drupke commented 1 year ago

OK, I fixed the radius=0 thing, which now defaults to a single-spaxel extraction.

WujiWang-astro commented 1 year ago

Hi Dave,

Thank you! I am not quite familiar with how GitHub coding works. I just checked the source code, it seems that
if radius ==0 then radius = 0.01 is still there. But this is a minor thing that is not why I see different flux levels of the order of 1e-4 (see below). More onto the spectra extraction part, maybe it gives more flexibility for q3dfit when extracting the spectrum when we can also specify the method of photutils.aperture.aperture_photometry, like method='center' or method='subpixel'. Again, this is not important now.

The thing I finally found that bothers me about the 1e-4 is in the read.py, line 259: the code tries to see if there is the keyword for the unit 'BUNIT' in the header. This may cause trouble but not necessarily give error message (like the case for me). For example, if the 'BUNIT' in the header is not correct, or the keyword is there but nothing (blank string) under the header['BUNIT'], the user cannot overwrite with 'fluxunit_in' from q3dift, but have to correct the header first. Maybe just add an output message there like 'q3dfit finds the unit xxx in the header and ignores the users input' such that it is easier for the user to know what is going on.

Wuji

drupke commented 1 year ago

Hi Wuji,

if radius ==0 then radius = 0.01 is still there.

Sorry, I should have clarified that this is pushed to the q3dfit-dev main branch, but not yet to the public q3dfit main. I want to do a little more testing before another release. If you can’t access the latest q3dfit-dev branch, let me know and I can give you access. But this is a minor thing that is not why I see different flux levels of the order of 1e-4 (see below). More onto the spectra extraction part, maybe it gives more flexibility for q3dfit when extracting the spectrum when we can also specify the method of photutils.aperture.aperture_photometry, like method='center' or method='subpixel'. Again, this is not important now.

The thing I finally found that bothers me about the 1e-4 is in the read.py, line 259: the code tries to see if there is the keyword for the unit 'BUNIT' in the header. This may cause trouble but not necessarily give error message (like the case for me). For example, if the 'BUNIT' in the header is not correct, or the keyword is there but nothing (blank string) under the header['BUNIT'], the user cannot overwrite with 'fluxunit_in' from q3dift, but have to correct the header first. Maybe just add an output message there like 'q3dfit finds the unit xxx in the header and ignores the users input' such that it is easier for the user to know what is going on.

OK, this is helpful. I’ll make some sort of fix for this so that BUNIT can be made optional. Will keep you posted.

Dave