Multi-Order Coverage maps represent regions in a spheric surface defined by tree-like structures with the aim of producing maps through different spatial resolutions. In astronomy they are used to procude lightweight version of surveys' coverage.
MOCs build the maps through Healpix. Healpix/MOC maps split the sky surface in diamond-like figures covering the same area each. Initially, at level zero, the number of diamonds is 12, covering the sphere in diamonds sized ~58 degree. At each following level, the diamonds are sub-divide in 4, recursively. Until the maximum number of levels is reached, each level doubles the map resolution, splitting the sphere firstly in 48 diamond, at the second level in 192; third, 768 and so on. The maximum level -- or order -- is 29, covering the sphere with 393.2 microarcsecond size diamonds.
In [1]:
# Let's handle units
from astropy import units as u
# Structure to map healpix' levels to their angular sizes
#
healpix_levels = {
0 : 58.63 * u.deg,
1 : 29.32 * u.deg,
2 : 14.66 * u.deg,
3 : 7.329 * u.deg,
4 : 3.665 * u.deg,
5 : 1.832 * u.deg,
6 : 54.97 * u.arcmin,
7 : 27.48 * u.arcmin,
8 : 13.74 * u.arcmin,
9 : 6.871 * u.arcmin,
10 : 3.435 * u.arcmin,
11 : 1.718 * u.arcmin,
12 : 51.53 * u.arcsec,
13 : 25.77 * u.arcsec,
14 : 12.88 * u.arcsec,
15 : 6.442 * u.arcsec,
16 : 3.221 * u.arcsec,
17 : 1.61 * u.arcsec,
18 : 805.2 * u.milliarcsecond,
19 : 402.6 * u.milliarcsecond,
20 : 201.3 * u.milliarcsecond,
21 : 100.6 * u.milliarcsecond,
22 : 50.32 * u.milliarcsecond,
23 : 25.16 * u.milliarcsecond,
24 : 12.58 * u.milliarcsecond,
25 : 6.291 * u.milliarcsecond,
26 : 3.145 * u.milliarcsecond,
27 : 1.573 * u.milliarcsecond,
28 : 786.3 * u.microarcsecond,
29 : 393.2 * u.microarcsecond
}
The libraries we can use to generate/manipulate Healpix/MOC maps are:
In [2]:
# as usual, matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
# Load Healpix
import healpy
# Erin Sheldon's healpix_util
import healpix_util as hu
# Thomas Boch's MOCpy
import mocpy
In [3]:
%ls
In [4]:
from astropy.io import fits
chandra = fits.open('Chandra_multiwavelength.fits')[1]
In [5]:
print chandra.columns
In [6]:
# we are interested here on columns 'RA','DEC' and 'RADEC_ERR'
_data = {'ra' : chandra.data['RA'] * u.degree,
'dec': chandra.data['DEC']* u.degree,
'pos_err' : chandra.data['RADEC_ERR']* u.arcsec}
In [7]:
from astropy.table import Table
_table = Table(_data)
In [8]:
import pandas as pd
df = _table.to_pandas()
del _table,_data
In [9]:
df.hist('pos_err',bins=100)
plt.show()
In [10]:
df.describe()
Out[10]:
In [11]:
# A function to find out which healpix level corresponds a given (typical) size of coverage
def size2level(size):
"""
Returns nearest Healpix level corresponding to a given diamond size
The 'nearest' Healpix level is here to be the nearest greater level,
right before the first level smaller than 'size'.
"""
assert size.unit
ko = None
for k,v in healpix_levels.iteritems():
if v < 2 * size:
break
ko = k
return ko
In [12]:
level = size2level(df.pos_err.median()* u.arcsec)
nside = 2**level
In [13]:
print "Typical (median) position error: \n{}".format(df.pos_err.median())
print "\nCorrespondig healpix level: {} \n\t and nsize value: {}".format(level,nside)
In [14]:
# Let's convert from ra,dec to theta,phi
# This function comes from mocpy
def ra2phi(ra):
"""
convert equatorial ra, dec in degrees
to polar theta, phi in radians
"""
import math
return math.radians(ra)
def dec2theta(dec):
"""
convert equatorial ra, dec in degrees
to polar theta, phi in radians
"""
import math
return math.pi/2 - math.radians(dec)
def radec2thetaphi(ra,dec):
_phi = ra2phi(ra)
_theta = dec2theta(dec)
return _theta,_phi
In [15]:
import healpy
def healpix_radec2pix(nside, ra, dec, nest=True):
_theta,_phi = radec2thetaphi(ra, dec)
return healpy.ang2pix(nside, _theta, _phi, nest=nest)
In [16]:
df['phi'] = df.ra.apply(ra2phi)
df['theta'] = df.dec.apply(dec2theta)
In [17]:
df.describe()
Out[17]:
In [18]:
hp_pix_eq = df.apply(lambda x:healpix_radec2pix(nside,x.ra,x.dec,nest=True), axis=1)
hp_pix_ang = df.apply(lambda x:healpy.ang2pix(nside,x.theta,x.phi,nest=True), axis=1)
import numpy
numpy.array_equal(hp_pix_ang,hp_pix_eq)
Out[18]:
In [19]:
hpix = hu.HealPix(scheme='nest',nside=nside)
hpix
Out[19]:
In [20]:
hu_pix = hpix.eq2pix(ra=df.ra,dec=df.dec)
numpy.array_equal(hu_pix,hp_pix_ang) and numpy.array_equal(hu_pix,hp_pix_eq)
Out[20]:
In [21]:
# Curiosity: which one is faster?
%timeit hpix.eq2pix(ra=df.ra,dec=df.dec)
%timeit df.apply(lambda x:healpix_radec2pix(nside,x.ra,x.dec,nest=True), axis=1)
In [22]:
# So...all results are equal \o/ and ES's is faster
# we can now go on and put it inside our DataFrame
df['hpix'] = hu_pix
In [23]:
df.describe()
Out[23]:
In [24]:
moc = mocpy.MOC()
moc.add_pix_list(level,df.hpix)
In [25]:
moc.plot()
In [26]:
moc.write('chandra_MOC_uniq.fits')
In [27]:
table = Table.from_pandas(df)
table.write('chandra_MOC_radec.fits',format='fits',overwrite=True)
In [28]:
del df,table,moc,chandra,hpix
%ls -lh
In [29]:
from astropy.io import fits
xmm = fits.open('XMM_multiwavelength_cat.fits')[1]
xmm.columns.names
Out[29]:
In [30]:
# we are interested here on columns 'RA','DEC' and 'RADEC_ERR'
_data = {'ra' : xmm.data['RA'] * u.degree,
'dec': xmm.data['DEC']* u.degree,
'pos_err' : xmm.data['RADEC_ERR']* u.arcsec}
df = Table(_data).to_pandas()
df.hist('pos_err',bins=100)
plt.show()
df.describe()
Out[30]:
In [31]:
level = size2level(df.pos_err.median()* u.arcsec)
nside = 2**level
print "Typical (median) position error: \n{}".format(df.pos_err.median())
print "\nCorrespondig healpix level: {} \n\t and nsize value: {}".format(level,nside)
In [32]:
hpix = hu.HealPix(scheme='nest',nside=nside)
hpix
Out[32]:
In [33]:
df['hpix'] = hpix.eq2pix(ra=df.ra,dec=df.dec)
df.describe()
Out[33]:
In [34]:
moc = mocpy.MOC()
moc.add_pix_list(level,df.hpix)
moc.plot()
In [35]:
moc.write('xmm_MOC_uniq.fits')
table = Table.from_pandas(df)
table.write('xmm_MOC_radec.fits',format='fits',overwrite=True)
In [36]:
%ls -lh
We need, to generate a coverage map, the following parameters:
Object::
given params:
_ra -> <vector>
_dec -> <vector>
_radius -> <vector>
computed params:
_healpix-level
_healpix-nside
output:
_moc-uniq <- <filename>
_moc-obj <- <filename>
In [37]:
def radec_2_moc(filename,ra_column,dec_column,radius_column=None,radius_value=None):
import healpix_util
import mocpy
import time
start_all = time.clock()
tbhdu = open_fits(filename)
table = radec_table(tbhdu,ra_column,dec_column,radius_column)
start_convert = time.clock()
if not radius_column:
if radius_value != None and radius_value > 0:
radius = radius_value
else:
from astropy import units
radius = 1 * units.arcsec
else:
radius = radius_mean(tbhdu,radius_column)
assert hasattr(radius,'unit')
level = size2level(radius)
nside = 2**level
hpix = healpix_util.HealPix('nest',nside)
table['hpix'] = hpix.eq2pix(table['ra'],table['dec'])
stop_convert = time.clock()
fileroot = '.'.join(filename.split('.')[:-1])
start_write_normal = time.clock()
fileout = '_'.join([fileroot,'MOC_position.fit'])
table.write(fileout,format='fits',overwrite=True)
stop_write_normal = time.clock()
start_write_moc = time.clock()
# fileout = '_'.join([fileroot,'MOC_uniq.fit'])
# moc = mocpy.MOC()
# moc.add_pix_list(level,table['hpix'])
# moc.write(fileout)
stop_write_moc = time.clock()
stop_all = time.clock()
_msg = "Time elapsed converting pixels: {}\n".format(stop_convert-start_convert)
_msg += "Time elapsed on writing the table: {}\n".format(stop_write_normal-start_write_normal)
_msg += "Time elapsed on writing MOC: {}\n".format(stop_write_moc-start_write_moc)
_msg += "Total time: {}\n".format(stop_all-start_all)
_msg += "Number of points: {}\n".format(len(table))
return _msg
def open_fits(filename,hdu=1):
from astropy.io import fits
from astropy.units import Quantity
_tab = fits.open(filename,ignore_missing_end=True)[hdu]
return _tab
def radec_table(tbhdu,ra_column,dec_column,radius_column=None):
from astropy.table import Table
from astropy import units
import numpy
_data = {'ra':tbhdu.data.field(ra_column) * units.deg,
'dec':tbhdu.data.field(dec_column) * units.deg,
'id':numpy.arange(tbhdu.header['NAXIS2'])}
if radius_column:
try:
_d = tbhdu.data.field(radius_column)
_data.update({'radius':_d})
except:
pass
return Table(_data)
def radius_mean(tbhdu,radius_column):
from astropy.units import Quantity
radius = None
if radius_column:
_radius = Quantity(tbhdu.data.field(radius_column), u.arcsec)
radius = _radius.mean()
assert radius
return radius
In [38]:
res = radec_2_moc('Chandra_multiwavelength.fits','RA','DEC','RADEC_ERR')
print res
In [39]:
res = radec_2_moc('XMM_multiwavelength_cat.fits','RA','DEC','RADEC_ERR')
print res
In [40]:
%ls -lh
In [41]:
def print_fits_columns(fitsfile,hdu=1):
from astropy.io import fits
hdul = fits.open(fitsfile,ignore_missing_end=True)
tbhdu = hdul[1]
print "Number of objects: {}\n".format(tbhdu.header['NAXIS2'])
print "{} columns:\n".format(fitsfile)
ncols = len(tbhdu.columns)
i = 0
for c in tbhdu.columns:
if i<=5:
print "\t{}; ".format(c.name)
else:
print "\t... ({} columns)".format(ncols-i)
break
i += 1
hdul.close()
In [42]:
print_fits_columns('photometry/hers/hers_catalogue_3sig250_no_extended.fits')
In [43]:
res = radec_2_moc('photometry/hers/hers_catalogue_3sig250_no_extended.fits','RA','DEC')
print res
In [44]:
print_fits_columns('photometry/galex/S82_gmsc_chbrandt.fit')
In [45]:
res = radec_2_moc('photometry/galex/S82_gmsc_chbrandt.fit','ra','dec','poserr')
print res
In [46]:
print_fits_columns('photometry/sdss/Stripe82_photo_chbrandt.fit')
In [47]:
res = radec_2_moc('photometry/sdss/Stripe82_photo_chbrandt.fit','ra','dec')
print res
In [48]:
print_fits_columns('photometry/shela/shela_stripe82_v1.3_cat.fits')
In [49]:
res = radec_2_moc('photometry/shela/shela_stripe82_v1.3_cat.fits','SDSS_RA','SDSS_DEC')
print res
In [50]:
print_fits_columns('photometry/spies/SpIES_ch1ch2_allaor_5s_bothchan_final.fits')
In [51]:
res = radec_2_moc('photometry/spies/SpIES_ch1ch2_allaor_5s_bothchan_final.fits','RA','DEC')
print res
In [52]:
print_fits_columns('photometry/unwise/brandt.fits')
In [53]:
res = radec_2_moc('photometry/unwise/brandt.fits','ra','dec')
print res
In [54]:
print_fits_columns('photometry/vla/first_14dec17.fits')
In [55]:
res = radec_2_moc('photometry/vla/first_14dec17.fits','RA','DEC')
print res
In [ ]: