Querying ALMA archive for M83 pointings and plotting them on a 2MASS image
from astroquery.alma import Alma
from astroquery.skyview import SkyView
import string
from astropy import units as u
import pylab as pl
import aplpy
Retrieve M83 2MASS K-band image:
m83_images = SkyView.get_images(position='M83', survey=['2MASS-K'], pixels=1500)
WARNING: AstropyDeprecationWarning: Since 0.4, config parameter 'astropy.utils.data.REMOTE_TIMEOUT' is deprecated. Use 'astropy.utils.data.conf.remote_timeout' instead. [astropy.config.configuration] WARNING:astropy:AstropyDeprecationWarning: Since 0.4, config parameter 'astropy.utils.data.REMOTE_TIMEOUT' is deprecated. Use 'astropy.utils.data.conf.remote_timeout' instead. FITS files must be read as binaries; error is likely.
Retrieve ALMA archive information including private data and non-science fields:
m83 = Alma.query_object('M83', public=False, science=False)
m83
Project code | Source name | RA | Dec | Band | Frequency resolution | Integration | Release date | Frequency support | Velocity resolution | Pol products | Observation date | PI name | PWV | Member ous id | Asdm uid | Project title | Project type | Scan intent |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
deg | deg | kHz | s | GHz | m / s | mm | ||||||||||||
string128 | string2048 | float64 | float64 | int32 | float64 | float64 | object | string21600 | float64 | string512 | object | string2048 | float32 | string256 | string256 | string2048 | string128 | string2048 |
2011.0.00772.S | M83 | 204.24697661723459 | -29.854173214826339 | 3 | 488.28125 | 40.034999999999997 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-03-27 07:12:59 | Hirota, Akihiko | 1.3226013 | uid://A002/X3216af/X31 | uid://A002/X3b3400/X90f | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24697661723459 | -29.854173214826339 | 3 | 488.28125 | 40.034999999999997 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-03-27 08:41:25 | Hirota, Akihiko | 1.2252344 | uid://A002/X3216af/X31 | uid://A002/X3b3400/Xaf3 | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24697661723459 | -29.854173214826339 | 3 | 488.28125 | 40.173000000000002 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-05-07 02:48:00 | Hirota, Akihiko | 1.8213283 | uid://A002/X3216af/X31 | uid://A002/X3fbe56/X607 | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24848880184393 | -29.864865672087014 | 3 | 488.28125 | 46.731000000000002 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-03-27 07:12:59 | Hirota, Akihiko | 1.3226013 | uid://A002/X3216af/X31 | uid://A002/X3b3400/X90f | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24848880184393 | -29.864865672087014 | 3 | 488.28125 | 46.75 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-03-27 08:41:25 | Hirota, Akihiko | 1.2252344 | uid://A002/X3216af/X31 | uid://A002/X3b3400/Xaf3 | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24848880184393 | -29.864865672087014 | 3 | 488.28125 | 46.909999999999997 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-05-07 02:48:00 | Hirota, Akihiko | 1.8213283 | uid://A002/X3216af/X31 | uid://A002/X3fbe56/X607 | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24977987250298 | -29.848448429708835 | 3 | 488.28125 | 40.030999999999999 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-03-27 07:12:59 | Hirota, Akihiko | 1.3226013 | uid://A002/X3216af/X31 | uid://A002/X3b3400/X90f | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24977987250298 | -29.848448429708835 | 3 | 488.28125 | 40.036000000000001 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-03-27 08:41:25 | Hirota, Akihiko | 1.2252344 | uid://A002/X3216af/X31 | uid://A002/X3b3400/Xaf3 | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
2011.0.00772.S | M83 | 204.24977987250298 | -29.848448429708835 | 3 | 488.28125 | 40.179000000000002 | 2013-09-28 | [100.63..101.57GHz,488.28kHz, XX YY] U [102.43..103.37GHz,488.28kHz, XX YY] U [112.74..113.68GHz,488.28kHz, XX YY] U [114.45..115.38GHz,488.28kHz, XX YY] | 1354.7258237091098 | XX YY | 2012-05-07 02:48:00 | Hirota, Akihiko | 1.8213283 | uid://A002/X3216af/X31 | uid://A002/X3fbe56/X607 | Giant Molecular Cloud Survey Toward bar and arm of the nearby Galaxy M83 | S | TARGET |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
2012.1.00762.S | m83 | 204.29994502155807 | -29.864462345910045 | 3 | 976.5625 | 10.129 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7724487100772 | XX YY | 2013-12-04 09:35:37 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X74fea5/X10d0 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.29994774826423 | -29.864447897690617 | 3 | 976.5625 | 10.151999999999999 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7823942972855 | XX YY | 2013-12-01 09:19:08 | Hirota, Akihiko | 0.54483771 | uid://A002/X5a9a13/X68b | uid://A002/X74d4df/X4b9 | Extended GMC survey in the nearby galaxy M83 | S | TARGET WVR | |
2012.1.00762.S | m83 | 204.2999565036323 | -29.87677639346035 | 3 | 976.5625 | 10.117000000000001 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7343188001446 | XX YY | 2013-12-17 08:47:59 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X75f169/X1348 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.30003398242312 | -29.833753778779368 | 3 | 976.5625 | 10.117000000000001 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7352830617078 | XX YY | 2013-12-17 10:44:34 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X75f169/X15f6 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.30007748596259 | -29.858283377392969 | 3 | 976.5625 | 10.119 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7789019763472 | XX YY | 2013-12-02 08:41:32 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X74deb6/X1304 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.31232521341892 | -29.815279582885864 | 3 | 976.5625 | 10.128 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7352830617078 | XX YY | 2013-12-17 10:44:34 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X75f169/X15f6 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.31233028180358 | -29.827583288376086 | 3 | 976.5625 | 10.118 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7732219137765 | XX YY | 2013-12-04 11:19:47 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X74fea5/X1558 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.31236333974948 | -29.839846902175914 | 3 | 976.5625 | 10.129 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7724487100772 | XX YY | 2013-12-04 09:35:37 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X74fea5/X10d0 | Extended GMC survey in the nearby galaxy M83 | S | TARGET | |
2012.1.00762.S | m83 | 204.31237002788583 | -29.852143686507858 | 3 | 976.5625 | 10.151999999999999 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7823942972855 | XX YY | 2013-12-01 09:19:08 | Hirota, Akihiko | 0.54483771 | uid://A002/X5a9a13/X68b | uid://A002/X74d4df/X4b9 | Extended GMC survey in the nearby galaxy M83 | S | TARGET WVR | |
2012.1.00762.S | m83 | 204.31237891503417 | -29.864472325574329 | 3 | 976.5625 | 10.117000000000001 | [99.86..101.85GHz,976.56kHz, XX YY] U [101.74..103.73GHz,976.56kHz, XX YY] U [112.09..114.08GHz,976.56kHz, XX YY] U [114.59..115.59GHz,488.28kHz, XX YY] | 2611.7343188001446 | XX YY | 2013-12-17 08:47:59 | Hirota, Akihiko | -- | uid://A002/X5a9a13/X68b | uid://A002/X75f169/X1348 | Extended GMC survey in the nearby galaxy M83 | S | TARGET |
Parse components of the ALMA data. Specifically, find the frequency support - the frequency range covered - and convert that into a central frequency for beam radius estimation.
def parse_frequency_support(frequency_support_str):
supports = frequency_support_str.split("U")
freq_ranges = [(float(sup.strip('[] ').split("..")[0]),
float(sup.strip('[] ').split("..")[1].split(',')[0].strip(string.letters)))
*u.Unit(sup.strip('[] ').split("..")[1].split(',')[0].strip(string.punctuation+string.digits))
for sup in supports]
return u.Quantity(freq_ranges)
def approximate_primary_beam_sizes(frequency_support_str):
freq_ranges = parse_frequency_support(frequency_support_str)
beam_sizes = [(1.22*fr.mean().to(u.m, u.spectral())/(12*u.m)).to(u.arcsec,
u.dimensionless_angles())
for fr in freq_ranges]
return u.Quantity(beam_sizes)
primary_beam_radii = [approximate_primary_beam_sizes(row['Frequency support']) for row in m83]
Compute primary beam parameters for the public and private components of the data for plotting below.
print "The bands used include: ",np.unique(m83['Band'])
The bands used include: Band ---- 3 6
private_circle_parameters = [(row['RA'],row['Dec'],np.mean(rad).to(u.deg).value)
for row,rad in zip(m83, primary_beam_radii)
if row['Release date']!='' and row['Band']==3]
public_circle_parameters = [(row['RA'],row['Dec'],np.mean(rad).to(u.deg).value)
for row,rad in zip(m83, primary_beam_radii)
if row['Release date']=='' and row['Band']==3]
unique_private_circle_parameters = np.array(list(set(private_circle_parameters)))
unique_public_circle_parameters = np.array(list(set(public_circle_parameters)))
print "BAND 3"
print "PUBLIC: Number of rows: {0}. Unique pointings: {1}".format(len(m83), len(unique_public_circle_parameters))
print "PRIVATE: Number of rows: {0}. Unique pointings: {1}".format(len(m83), len(unique_private_circle_parameters))
private_circle_parameters_band6 = [(row['RA'],row['Dec'],np.mean(rad).to(u.deg).value)
for row,rad in zip(m83, primary_beam_radii)
if row['Release date']!='' and row['Band']==6]
public_circle_parameters_band6 = [(row['RA'],row['Dec'],np.mean(rad).to(u.deg).value)
for row,rad in zip(m83, primary_beam_radii)
if row['Release date']=='' and row['Band']==6]
BAND 3 PUBLIC: Number of rows: 2720. Unique pointings: 551 PRIVATE: Number of rows: 2720. Unique pointings: 650
Show all of the private observation pointings that have been acquired
fig = aplpy.FITSFigure(m83_images[0])
fig.show_grayscale(stretch='arcsinh')
fig.show_circles(unique_private_circle_parameters[:,0],
unique_private_circle_parameters[:,1],
unique_private_circle_parameters[:,2],
color='r', alpha=0.2)
/Users/adam/repos/aplpy/aplpy/labels.py:432: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if self.coord == x or self.axis.apl_tick_positions_world[ipos] > 0: INFO:astropy:Auto-setting vmin to 5.516e+02 INFO:astropy:Auto-setting vmax to 5.790e+02
INFO: Auto-setting vmin to 5.516e+02 [aplpy.core] INFO: Auto-setting vmax to 5.790e+02 [aplpy.core]
In principle, all of the pointings shown below should be downloadable from the archive:
fig = aplpy.FITSFigure(m83_images[0])
fig.show_grayscale(stretch='arcsinh')
fig.show_circles(unique_public_circle_parameters[:,0],
unique_public_circle_parameters[:,1],
unique_public_circle_parameters[:,2],
color='b', alpha=0.2)
INFO:astropy:Auto-setting vmin to 5.516e+02 INFO:astropy:Auto-setting vmax to 5.790e+02
INFO: Auto-setting vmin to 5.516e+02 [aplpy.core] INFO: Auto-setting vmax to 5.790e+02 [aplpy.core]
Use pyregion to write the observed regions to disk. Pyregion has a very awkward API; there is (in principle) work in progress to improve that situation but for now one must do all this extra work.
import pyregion
from pyregion.parser_helper import Shape
prv_regions = pyregion.ShapeList([Shape('circle',[x,y,r]) for x,y,r in private_circle_parameters])
pub_regions = pyregion.ShapeList([Shape('circle',[x,y,r]) for x,y,r in public_circle_parameters])
for r,(x,y,c) in zip(prv_regions+pub_regions,
np.vstack([private_circle_parameters,
public_circle_parameters])):
r.coord_format = 'fk5'
r.coord_list = [x,y,c]
r.attr = ([], {'color': 'green', 'dash': '0 ', 'dashlist': '8 3 ', 'delete': '1 ', 'edit': '1 ',
'fixed': '0 ', 'font': '"helvetica 10 normal roman"', 'highlite': '1 ',
'include': '1 ', 'move': '1 ', 'select': '1 ', 'source': '1', 'text': '',
'width': '1 '})
prv_regions.write('M83_observed_regions_private_March2015.reg')
pub_regions.write('M83_observed_regions_public_March2015.reg')
from astropy.io import fits
prv_mask = fits.PrimaryHDU(prv_regions.get_mask(m83_images[0][0]).astype('int'),
header=m83_images[0][0].header)
pub_mask = fits.PrimaryHDU(pub_regions.get_mask(m83_images[0][0]).astype('int'),
header=m83_images[0][0].header)
WARNING: AstropyDeprecationWarning: The ascard function is deprecated and may be removed in a future version. Use the `.cards` attribute instead. [astropy.utils.decorators] WARNING:astropy:AstropyDeprecationWarning: The ascard function is deprecated and may be removed in a future version. Use the `.cards` attribute instead. WARNING: AstropyDeprecationWarning: The CardList class has been deprecated; all its former functionality has been subsumed by the Header class, so CardList objects should not be directly created. See the PyFITS 3.1.0 CHANGELOG for more details. [astropy.io.fits.card] WARNING:astropy:AstropyDeprecationWarning: The CardList class has been deprecated; all its former functionality has been subsumed by the Header class, so CardList objects should not be directly created. See the PyFITS 3.1.0 CHANGELOG for more details. WARNING: AstropyDeprecationWarning: The ascard function is deprecated and may be removed in a future version. Use the `.cards` attribute instead. [astropy.utils.decorators] WARNING:astropy:AstropyDeprecationWarning: The ascard function is deprecated and may be removed in a future version. Use the `.cards` attribute instead.
pub_mask.writeto('public_m83_almaobs_mask.fits', clobber=True)
fig = aplpy.FITSFigure(m83_images[0])
fig.show_grayscale(stretch='arcsinh')
fig.show_contour(prv_mask, levels=[0.5,1], colors=['r','r'])
fig.show_contour(pub_mask, levels=[0.5,1], colors=['b','b'])
INFO:astropy:Auto-setting vmin to 5.516e+02 INFO:astropy:Auto-setting vmax to 5.790e+02
INFO: Auto-setting vmin to 5.516e+02 [aplpy.core] INFO: Auto-setting vmax to 5.790e+02 [aplpy.core]
Now we create a 'hit mask' showing the relative depth of each observed field in each band
hit_mask_band3_public = np.zeros_like(m83_images[0][0].data)
hit_mask_band3_private = np.zeros_like(m83_images[0][0].data)
hit_mask_band6_public = np.zeros_like(m83_images[0][0].data)
hit_mask_band6_private = np.zeros_like(m83_images[0][0].data)
from astropy import wcs
mywcs = wcs.WCS(m83_images[0][0].header)
for row,rad in zip(m83, primary_beam_radii):
shape = Shape('circle', (row['RA'], row['Dec'],np.mean(rad).to(u.deg).value))
shape.coord_format = 'fk5'
shape.coord_list = (row['RA'], row['Dec'],np.mean(rad).to(u.deg).value)
shape.attr = ([], {'color': 'green', 'dash': '0 ', 'dashlist': '8 3 ', 'delete': '1 ', 'edit': '1 ',
'fixed': '0 ', 'font': '"helvetica 10 normal roman"', 'highlite': '1 ',
'include': '1 ', 'move': '1 ', 'select': '1 ', 'source': '1', 'text': '',
'width': '1 '})
if row['Release date']=='' and row['Band']==3:
(xlo,xhi,ylo,yhi),mask = pyregion_subset(shape, hit_mask_band3_private, mywcs)
hit_mask_band3_private[ylo:yhi,xlo:xhi] += row['Integration']*mask
elif row['Release date'] and row['Band']==3:
(xlo,xhi,ylo,yhi),mask = pyregion_subset(shape, hit_mask_band3_public, mywcs)
hit_mask_band3_public[ylo:yhi,xlo:xhi] += row['Integration']*mask
elif row['Release date'] and row['Band']==6:
(xlo,xhi,ylo,yhi),mask = pyregion_subset(shape, hit_mask_band6_public, mywcs)
hit_mask_band6_public[ylo:yhi,xlo:xhi] += row['Integration']*mask
elif row['Release date']=='' and row['Band']==6:
(xlo,xhi,ylo,yhi),mask = pyregion_subset(shape, hit_mask_band6_private, mywcs)
hit_mask_band6_private[ylo:yhi,xlo:xhi] += row['Integration']*mask
fig = aplpy.FITSFigure(m83_images[0])
fig.show_grayscale(stretch='arcsinh')
fig.show_contour(fits.PrimaryHDU(data=hit_mask_band3_public, header=m83_images[0][0].header),
levels=np.logspace(0,5,base=2, num=6), colors=['r']*6)
fig.show_contour(fits.PrimaryHDU(data=hit_mask_band3_private, header=m83_images[0][0].header),
levels=np.logspace(0,5,base=2, num=6), colors=['y']*6)
fig.show_contour(fits.PrimaryHDU(data=hit_mask_band6_public, header=m83_images[0][0].header),
levels=np.logspace(0,5,base=2, num=6), colors=['c']*6)
fig.show_contour(fits.PrimaryHDU(data=hit_mask_band6_private, header=m83_images[0][0].header),
levels=np.logspace(0,5,base=2, num=6), colors=['b']*6)
INFO:astropy:Auto-setting vmin to 5.516e+02 INFO:astropy:Auto-setting vmax to 5.790e+02
INFO: Auto-setting vmin to 5.516e+02 [aplpy.core] INFO: Auto-setting vmax to 5.790e+02 [aplpy.core]
from astropy import wcs
import pyregion
from astropy import log
def pyregion_subset(region, data, mywcs):
"""
Return a subset of an image (`data`) given a region.
"""
shapelist = pyregion.ShapeList([region])
if shapelist[0].coord_format not in ('physical','image'):
# Requires astropy >0.4...
# pixel_regions = shapelist.as_imagecoord(self.wcs.celestial.to_header())
# convert the regions to image (pixel) coordinates
celhdr = mywcs.sub([wcs.WCSSUB_CELESTIAL]).to_header()
pixel_regions = shapelist.as_imagecoord(celhdr)
else:
# For this to work, we'd need to change the reference pixel after cropping.
# Alternatively, we can just make the full-sized mask... todo....
raise NotImplementedError("Can't use non-celestial coordinates with regions.")
pixel_regions = shapelist
# This is a hack to use mpl to determine the outer bounds of the regions
# (but it's a legit hack - pyregion needs a major internal refactor
# before we can approach this any other way, I think -AG)
mpl_objs = pixel_regions.get_mpl_patches_texts()[0]
# Find the minimal enclosing box containing all of the regions
# (this will speed up the mask creation below)
extent = mpl_objs[0].get_extents()
xlo, ylo = extent.min
xhi, yhi = extent.max
all_extents = [obj.get_extents() for obj in mpl_objs]
for ext in all_extents:
xlo = xlo if xlo < ext.min[0] else ext.min[0]
ylo = ylo if ylo < ext.min[1] else ext.min[1]
xhi = xhi if xhi > ext.max[0] else ext.max[0]
yhi = yhi if yhi > ext.max[1] else ext.max[1]
log.debug("Region boundaries: ")
log.debug("xlo={xlo}, ylo={ylo}, xhi={xhi}, yhi={yhi}".format(xlo=xlo,
ylo=ylo,
xhi=xhi,
yhi=yhi))
subwcs = mywcs[ylo:yhi, xlo:xhi]
subhdr = subwcs.sub([wcs.WCSSUB_CELESTIAL]).to_header()
subdata = data[ylo:yhi, xlo:xhi]
mask = shapelist.get_mask(header=subhdr,
shape=subdata.shape)
log.debug("Shapes: data={0}, subdata={2}, mask={1}".format(data.shape, mask.shape, subdata.shape))
return (xlo,xhi,ylo,yhi),mask