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 https://github.com/dendrograms/astrodendro/issues/139#issuecomment-100227139.

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:schilke@ph1.uni-koeln.de

[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 https://github.com/dendrograms/astrodendro/issues/139#issuecomment-100230819.

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:schilke@ph1.uni-koeln.de

[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 > https://github.com/dendrograms/astrodendro/issues/139#issuecomment-100236608. ## 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:schilke@ph1.uni-koeln.de

[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| , 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 https://github.com/dendrograms/astrodendro/issues/139#issuecomment-100240284.

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:schilke@ph1.uni-koeln.de

[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 https://github.com/dendrograms/astrodendro/issues/139#issuecomment-100243837.

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:schilke@ph1.uni-koeln.de