dendrograms / astrodendro

Generate a dendrogram from a dataset
https://dendrograms.readthedocs.io/
Other
37 stars 36 forks source link

PPVStatistic.flux wrongly ignores velocity units #107

Open tomr-stargazer opened 10 years ago

tomr-stargazer commented 10 years ago

Hi, I was working on getting integrated fluxes of PPV structures in a datacube using ppv_catalog. The unit handling for summing up the flux when the input units are (u.Jy/u.beam) seems to correctly handle integrating over 2d spatial units to produce u.Jy output, but the integration over the velocity axis should also contribute a velocity units term (usually I use u.km/u.s) to the computed flux (producing u.Jy * u.km / u.s) in order to be a physically meaningful sum over all 3 axes.

Currently I am doing a clumsy workaround that looks like this:

catalog = astrodendro.ppv_catalog(d, metadata)
flux_kms = catalog['flux'] * metadata['velocity_scale']

but it seems to me that the PPVStatistic and compute_flux code should be updated to correctly deal with this use case. Do you agree?

astrofrog commented 10 years ago

I have to think about this - if you have a cube where all pixels are in Jy/pixel then the total sum over all axes is going to be in Jy too. It will only be in Jy * km/s if you multiply each pixel value with dv (in km/s). This is similar to the fact that in a 2-d image, summing pixels in Jy/pixel does not give Jy * sr. So I think the correct behavior may actually be ok. On the other hand, if the original units are Jy/(km/s) then we should multiply by dv for each pixel to get Jy (like we do when we have Jy/beam and multiply by the beam area).

At the end of the day, compute_flux should always return Jy, and should do whatever is needed to be able to do that.

ChrisBeaumont commented 10 years ago

I had to look at http://www.astro.ufl.edu/~jt/alma/FLA_CDE_intro.pdf to refresh the units in my mind:

Jy is energy / area / time / frequency K = energy / area / time /frequency / solid angle

In my experience the most common unit for radio data is K or Jy/beam (both equivalent units). I think in these cases the integrated unit should thus be equivalent to K sr km/s, or Jy km/s

astrofrog commented 10 years ago

I'm managing to confuse myself. Say you have a ppv cube where each value is in Jy/pixel, and you want to make a pp image in Jy/pixel, what would you then do?

keflavich commented 10 years ago

Average it along the spectral axis; Jy is already per Hz. Or take the peak. Not the sum or integral unless you wanted Jy Hz / K km/s

ChrisBeaumont commented 10 years ago

I agree with Adam. Although K km/s is a common integrated flux measure, so there are use cases for summing over the third axis, and slapping on extra spectral unit

astrofrog commented 10 years ago

I guess what I would like compute_flux to do is ultimately return a value in Jy for all structures, since that can then be easily converted to mass, which is what we usually care about. To me it doesn't make sense to give the flux of a structure in e.g. K km/s. So I feel that we should change how the flux is combined rather than changing the units. Does that make some sense?

tomr-stargazer commented 10 years ago

@astrofrog @keflavich The use-case I have is to try and calculate the luminosity of various features, which requires that the datacube has been integrated in both solid angle and velocity to produce something that has K * km/s, as in the following equation: image (citation @low-sky et al 2008 -- http://adsabs.harvard.edu/abs/2008ApJ...679.1338R)

so @astrofrog it is not the case that Jy can be easily converted to mass unless you are using a different approach than the one I just cited.

keflavich commented 10 years ago

Technically I think K km/s (or Jy Hz) is flux while Jy or K is flux_density (i.e., flux per Hz). Chris, feel free to contradict me if you think I'm mistaken. So it may make sense to rename the function compute_flux_density and have the return unit be Jy, since as you say that is most readily convertible to mass.

keflavich commented 10 years ago

It really does make sense to convert K km/s to column/mass for a spectral line, since the flux density at any given spectral pixel is not representative of the line of sight. For dust, on the other hand, you really want flux density in Jy because that does represent the mass.

ChrisBeaumont commented 10 years ago

The equation above is the example I have in mind where K km/s comes up. Adam might be right about terms, which always confuse me (specific flux, specific intensity, flux, intensity, flux density, ...!). From an API perspective, I consider all these terms "broken" in the sense that keeping them all conceptually separated is harder than it needs to be. (Same thing with all the "temperature" terms, but that's a rant for another time)

Ultimately, both integrating and averaging over the third axis are reasonable operations, and the unit behavior of both is well defined. Can we rename 'flux' to something less ambiguous?

tomr-stargazer commented 10 years ago

@keflavich My understanding is also that Jy is "flux density" as noted in https://en.wikipedia.org/wiki/Jansky : image "It is important to understand the meaning of the per hertz component of the jansky unit. When measuring broadband continuum emissions, where the energy is roughly evenly distributed across the detector bandwidth, the detected signal will increase in proportion to the bandwidth of the detector (as opposed to signals with bandwidth narrower than the detector bandpass). To calculate the flux density in janskys, the total power detected (in watts) is divided by the receiver collecting area (in square meters), and then divided by the detector bandwidth (in hertz)."

I would agree with the idea of renaming compute_flux to compute_flux_density if a function that returns only Jy is still needed in a datacube like this. @ChrisBeaumont's comment about the ensuing API headache (perhaps we are already suffering from this) is worth considering.

low-sky commented 10 years ago

Right! This is the difference between continuum measures (for which Jy are appropriate) and spectral line measures (for which Jy Hz are appropriate). If you have a PPV data set, you are, by construction, in a spectral line regime.

Sorry to be terse; I have much more to say but have to go and be active in front of people.

e.

On Mar 27, 2014, at 9:09 AM, Adam Ginsburg notifications@github.com wrote:

It really does make sense to convert K km/s to column/mass for a spectral line, since the flux density at any given spectral pixel is not representative of the line of sight. For dust, on the other hand, you really want flux density in Jy because that does represent the mass.

— Reply to this email directly or view it on GitHub.

Erik Rosolowsky Asst. Professor • Dept. of Physics University of Alberta • +1-780-492-9272

astrofrog commented 10 years ago

Right, what's there currently should definitely be renamed to compute_flux_density. I was being dense on the Jy km/s or Jy Hz issue (which is a flux). Now the question is, do we then return a flux density if we are using a 2-d image, and a flux if we are using a 3-d image?

astrofrog commented 10 years ago

One option is to then apply compute_flux_density to all data, but only along the spatial dimensions. So for a 2-d image, we return a scalar, and for 3-d cubes we return an array of Jy values for each channel. Then in the 3-d case we have a routine that will sum up the 1-d array of Jy values along the spectral dimension to give a flux. Sounds reasonable?

astrofrog commented 10 years ago

Also, there is no problem in including a 1-d flux_density in the catalog table for PPV cubes - astropy's Table class supports vector columns, so we could definitely do that.

astrofrog commented 10 years ago

Also, in cases where flux_density is then a 1-d array, we simply also provide the integrated flux in the table. How does that sound?

keflavich commented 10 years ago

:+1:

tomr-stargazer commented 10 years ago

@astrofrog This approach is good! (where the flux is computed through the intermediate compute_flux_density step.)

Making any practical use of the vector of flux densities might be awkward in PPV data, since dendrogram structures usually have different spatial extent / number of pixels in each velocity channel (thus you are integrating over a different solid angle on each channel). But you can definitely sum them up.