Possible bug in subcube_from_ds9region?

[aakepley]

I used the subcube_from_ds9region code as a base to extract from regions from 2-d images. I noticed that the get_mpl_patches_texts method has a origin argument. I found I had to set this to origin=0 to get the appropriate pixels extracted, i.e.,

mpl_objs = region_image.get_mpl_patches_texts(origin=0)[0]

You may be taking care of this in the subcube method where you actually extract the data, but I'm still trying to parse exactly what's going on there. In that case, just ignore this issue.

[keflavich]

We call get_mask instead of get_mpl_patches_texts in spectral-cube: https://github.com/radio-astro-tools/spectral-cube/blob/master/spectral_cube/spectral_cube.py#L1804

The origin argument is either 1 or 0, depending on whether you are using a 0-indexed language (python, C) or a 1-indexed language (fortran, FITS).

Digging through the pyregion code, I actually can't tell if it's doing the right thing or not. origin=1, which is actually Wrong, is hard-coded in various places, which may be OK if it's done self-consistently, but I am worried now. I believe I have done tests where I compared the pyregion mask to the image directly in ds9 and it came out as expected, but it is worth double-checking.

p.s. keep an eye on https://github.com/astropy/regions, which will replace pyregion eventually

[aakepley]

Thanks for the quick response Adam! I'll definitely keep an eye on astropy/regions -- looks like a cool package.

For what it's worth, I did check my spectra by comparing spectra extracted manually via casaviewer and one extracted via spectral_cube and they did match. (Yay!)

I've included my 2d test code below in case it's useful to someone.

import numpy as np
import pyregion
from astropy.io import fits

# open image file and get header
fits_image = 'pacs70.reproject.fits'
image = fits.open(fits_image)
hdr = image[0].header
data = image[0].data

# open region list
region_file='origin_test.reg'
region_list =  pyregion.open(region_file)
region = region_list[0]

region_image = pyregion.ShapeList([region]).as_imagecoord(hdr)
mpl_objs = region_image.get_mpl_patches_texts(origin=0)[0]
extent = mpl_objs[0].get_extents()

xlo, ylo = extent.min
xhi, yhi = extent.max

data[ylo:yhi,xlo:xhi]