Unit of v_cen not given properly - no first moment

[PeterSchilke]

Hi,

I am looking at a PPV datacube, and the unit of v_cen is not given properly, just in pixels. But I also don't think this is the first moment (the intensity weighted velocity centroid), but just the header information, while v_rms is the proper second moment.

Peter

[keflavich]

Hi @PeterSchilke - are you looking at a catalog produced from the dendrogram? (How) have you specified the metadata?

[keflavich]

here is an example of metadata specification for a PPV cube: https://github.com/keflavich/APEX_CMZ_H2CO/blob/master/analysis/dendro_temperature.py#L60

[PeterSchilke]

I have specified the metadata - but they are not used, here the
definition of v_cen in analysis.py. The units are fine for v_rms, I
should have said that.

@Property
def v_cen(self):
"""
The mean velocity of the structure (where the velocity axis can be
specified by the vaxis metadata parameter, which defaults to 0
following the Numpy convention - the third axis in the FITS
convention).
"""
p = self._world_pos()
return p[self.vaxis]

On 05/08/2015 03:12 PM, Adam Ginsburg wrote:

Hi @PeterSchilke https://github.com/PeterSchilke - are you looking
at a catalog produced from the dendrogram? (How) have you specified
the metadata
http://www.dendrograms.org/en/latest/catalog.html?highlight=metadata#required-metadata?

—
Reply to this email directly or view it on GitHub
#139 (comment).

Prof. Dr. Peter Schilke
I. Physikalisches Institut der Universitaet zu Koeln
Zuelpicher Str. 77 D-50937 Koeln Germany
Tel. +49-221-470-1935 FAX +49-221-470-5162
mailto:[email protected]

[keflavich]

In principle, at least, _world_pos should be in velocity units (or frequency, or whatever the world units are). Since that's not happening, perhaps the WCS has not been properly attached to the cube. See https://github.com/dendrograms/astrodendro/blob/master/astrodendro/analysis.py#L294, especially the note:

            # TODO: set units correctly following WCS
            # We use origin=0 since the indices come from Numpy indexing
            return self.wcs.all_pix2world([xyz], 0).ravel()[::-1]

EDIT: this note implies that you might be getting the right values out of v_cen, i.e. in velocity, but with the wrong units.

[PeterSchilke]

This may well be - the file is read from fits, and if the WCS is not
attached automatically, I haven't done anything.

Anyway, the more interesting point for me would be to know if the first
moment can be calculated.

On 05/08/2015 03:21 PM, Adam Ginsburg wrote:

In principle, at least, |_world_pos| should be in velocity units (or
frequency, or whatever the world units are). Since that's not
happening, perhaps the WCS has not been properly attached to the cube.
See
https://github.com/dendrograms/astrodendro/blob/master/astrodendro/analysis.py#L294,
especially the note:

| # TODO: set units correctly following WCS
# We use origin=0 since the indices come from Numpy indexing
return self.wcs.all_pix2world([xyz], 0).ravel()[::-1]
|

—
Reply to this email directly or view it on GitHub
#139 (comment).

Prof. Dr. Peter Schilke
I. Physikalisches Institut der Universitaet zu Koeln
Zuelpicher Str. 77 D-50937 Koeln Germany
Tel. +49-221-470-1935 FAX +49-221-470-5162
mailto:[email protected]

[keflavich]

@PeterSchilke try something like this, if you've opened a FITS file:

from astropy.io import fits
from astropy import wcs
fitsfile = fits.open('filename.fits')

dend = Dendrogram.compute(fitsfile[0].data)

mywcs = wcs.WCS(fitsfile[0].header)
metadata = {}
metadata['wcs'] = mywcs

catalog = ppv_catalog(dend, metadata)

This isn't a complete example, but should give you a general idea of how to incorporate the WCS, which should then result in a correctly computed 1st moment

[PeterSchilke]

I do that, but get an error message:

In [39]: stat = PPVStatistic(d.trunk[0], metadata=metadata)

In [40]: stat.v_cen

ValueError Traceback (most recent call last)
in ()
----> 1 stat.v_cen

/usr/lib/python2.7/site-packages/astrodendro/analysis.pyc in v_cen(self)
399 following the Numpy convention - the third axis in the
FITS convention).
400 """
--> 401 p = self._world_pos()
402 return p[self.vaxis]
403

/usr/lib/python2.7/site-packages/astrodendro/analysis.pyc in
_world_pos(self)
263 # TODO: set units correctly following WCS
264 # We use origin=0 since the indices come from Numpy
indexing
--> 265 return self.wcs.all_pix2world([xyz], 0).ravel()[::-1]
266
267 @abc.abstractproperty

/usr/lib64/python2.7/site-packages/astropy/wcs/wcs.pyc in
all_pix2world(self, _args, *_kwargs)
1097 def all_pix2world(self, _args, *_kwargs):
1098 return self._array_converter(
-> 1099 self._all_pix2world, 'output', _args, *_kwargs)
1100 all_pix2world.doc = """
1101 Transforms pixel coordinates to world coordinates.

/usr/lib64/python2.7/site-packages/astropy/wcs/wcs.pyc in
_array_converter(self, func, sky, _args, *_kwargs)
1073 if self.naxis == 1 and len(xy.shape) == 1:
1074 return _return_list_of_arrays([xy], origin)
-> 1075 return _return_single_array(xy, origin)
1076
1077 elif len(args) == self.naxis + 1:

/usr/lib64/python2.7/site-packages/astropy/wcs/wcs.pyc in
_return_single_array(xy, origin)
1054 raise ValueError(
1055 "When providing two arguments, the array
must be "
-> 1056 "of shape (N, {0})".format(self.naxis))
1057 if ra_dec_order and sky == 'input':
1058 xy = self._denormalize_sky(xy)

ValueError: When providing two arguments, the array must be of shape (N, 4)

InOn 05/08/2015 03:38 PM, Adam Ginsburg wrote:

@PeterSchilke https://github.com/PeterSchilke try something like
this, if you've opened a FITS file:

|from astropy.io import fits
from astropy import wcs
fitsfile = fits.open('filename.fits')

dend = Dendrogram.compute(fitsfile[0].data)

mywcs = wcs.WCS(fitsfile[0].header)
metadata = {}
metadata['wcs'] = mywcs

catalog = ppv_catalog(dend, metadata)
|

This isn't a complete example, but should give you a general idea of
how to incorporate the WCS

—
Reply to this email directly or view it on GitHub
#139 (comment).

Prof. Dr. Peter Schilke
I. Physikalisches Institut der Universitaet zu Koeln
Zuelpicher Str. 77 D-50937 Koeln Germany
Tel. +49-221-470-1935 FAX +49-221-470-5162
mailto:[email protected]

[keflavich]

@PeterSchilke Ah, this makes sense! You have a 4D cube, where one dimension is the (degenerate) STOKES axis. There are two solutions you can try:

1. Instead of using astropy.io.fits to read the file, use spectral_cube, which will appropriately ditch the unnecessary axis
2. Run this command to get a more appropriate WCS with the STOKES axis dropped:

mywcs = wcs.WCS(fitsfile[0].header).sub([wcs.WCSSUB_CELESTIAL, wcs.WCSSUB_SPECTRAL])

[PeterSchilke]

I tried solution 2. No error message, but now the values have no units
at all.

In [77]: stat.y_cen
Out[77]: 0.78082086968815312

In [78]: stat.v_cen
Out[78]: -8570.7494935340383

In [79]: stat.x_cen
Out[79]: 351.23298147473946

WCSAXES = 3 / Number of coordinate
axes
CRPIX1 = 350.718126006 / Pixel coordinate of reference
point
CRPIX2 = 93.5508411697 / Pixel coordinate of reference
point
CRPIX3 = 186.0 / Pixel coordinate of reference
point
CDELT1 = -0.00397682344962 / [deg] Coordinate increment at reference
point
CDELT2 = 0.00397682344962 / [deg] Coordinate increment at reference
point
CDELT3 = 100.0 / [m/s] Coordinate increment at reference
point
CUNIT1 = 'deg' / Units of coordinate increment and
value
CUNIT2 = 'deg' / Units of coordinate increment and
value
CUNIT3 = 'm/s' / Units of coordinate increment and
value
CTYPE1 = 'GLON-SFL' / galactic longitude, Sanson-Flamsteed
projection
CTYPE2 = 'GLAT-SFL' / galactic latitude, Sanson-Flamsteed
projection
CTYPE3 = 'VOPT' / Optical velocity
(linear)
CRVAL1 = 351.180743055 / [deg] Coordinate value at reference
point
CRVAL2 = 0.7277389228 / [deg] Coordinate value at reference
point
CRVAL3 = 0.0 / [m/s] Coordinate value at reference
point
PV1_0 = 1.0 / [m/s] Coordinate value at reference
point
PV1_1 = 0.0 / [deg] Native longitude of the reference
point
PV1_2 = 0.7277389228 / [deg] Native latitude of the reference
point
LONPOLE = 0.0 / [deg] Native longitude of celestial
pole
LATPOLE = 90.0 / [deg] Native latitude of celestial
pole
RESTFRQ = 220398677000.0 / [Hz] Line rest
frequency
SPECSYS = 'LSRK' / Reference frame of spectral
coordinates
MJD-OBS = 51544.5 / [d] MJD of observation matching
DATE-OBS
DATE-OBS= '2000-01-01T12:00:00.0' / ISO-8601 observation date matching
MJD-OBS

On 05/08/2015 03:59 PM, Adam Ginsburg wrote:

@PeterSchilke https://github.com/PeterSchilke Ah, this makes sense!
You have a 4D cube, where one dimension is the (degenerate) STOKES
axis. There are two solutions you can try:

1. Instead of using |astropy.io.fits| to read the file, use
   |spectral_cube| <spectral-cube.rtfd.org>, which will appropriately
   ditch the unnecessary axis
2. Run this command to get a more appropriate WCS with the STOKES axis
   dropped:

|mywcs = wcs.WCS(fitsfile[0].header).sub([wcs.WCSSUB_CELESTIAL, wcs.WCSSUB_SPECTRAL])
|

—
Reply to this email directly or view it on GitHub
#139 (comment).

Prof. Dr. Peter Schilke
I. Physikalisches Institut der Universitaet zu Koeln
Zuelpicher Str. 77 D-50937 Koeln Germany
Tel. +49-221-470-1935 FAX +49-221-470-5162
mailto:[email protected]

[keflavich]

@PeterSchilke - yes, units for these are not implemented yet, but it looks like those values are in m/s

[PeterSchilke]

Thanks!

On 05/08/2015 04:11 PM, Adam Ginsburg wrote:

@PeterSchilke https://github.com/PeterSchilke - yes, units for these
are not implemented yet, but it looks like those values are in m/s

—
Reply to this email directly or view it on GitHub
#139 (comment).

Prof. Dr. Peter Schilke
I. Physikalisches Institut der Universitaet zu Koeln
Zuelpicher Str. 77 D-50937 Koeln Germany
Tel. +49-221-470-1935 FAX +49-221-470-5162
mailto:[email protected]