<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Processing MOC maps
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.
Step4:  The libraries we can use to generate/manipulate Healpix/MOC maps are
Step5:  Let's do the same with healpix_util now
Step6:  MOCpy for visualizing and writing the maps
Step7:  MOC catalogs now for LaMassa's XMM
Step8:  Automating MOC generation from fits files
We need, to generate a coverage map, the following parameters
<ASSISTANT_TASK:>
Python Code:
# 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
}
Explanation: Processing MOC maps
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.
End of explanation
# 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
%ls
from astropy.io import fits
chandra = fits.open('Chandra_multiwavelength.fits')[1]
print chandra.columns
# 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}
from astropy.table import Table
_table = Table(_data)
import pandas as pd
df = _table.to_pandas()
del _table,_data
df.hist('pos_err',bins=100)
plt.show()
df.describe()
# 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
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)
# 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
import healpy
def healpix_radec2pix(nside, ra, dec, nest=True):
    _theta,_phi = radec2thetaphi(ra, dec)
    return healpy.ang2pix(nside, _theta, _phi, nest=nest)
df['phi'] = df.ra.apply(ra2phi)
df['theta'] = df.dec.apply(dec2theta)
df.describe()
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)
Explanation: The libraries we can use to generate/manipulate Healpix/MOC maps are:
* healpy
* mocpy
* healpix_util
End of explanation
hpix = hu.HealPix(scheme='nest',nside=nside)
hpix
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)
# 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)
# 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
df.describe()
Explanation: Let's do the same with healpix_util now
End of explanation
moc = mocpy.MOC()
moc.add_pix_list(level,df.hpix)
moc.plot()
moc.write('chandra_MOC_uniq.fits')
table = Table.from_pandas(df)
table.write('chandra_MOC_radec.fits',format='fits',overwrite=True)
del df,table,moc,chandra,hpix
%ls -lh
Explanation: MOCpy for visualizing and writing the maps
End of explanation
from astropy.io import fits
xmm = fits.open('XMM_multiwavelength_cat.fits')[1]
xmm.columns.names
# 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()
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)
hpix = hu.HealPix(scheme='nest',nside=nside)
hpix
df['hpix'] = hpix.eq2pix(ra=df.ra,dec=df.dec)
df.describe()
moc = mocpy.MOC()
moc.add_pix_list(level,df.hpix)
moc.plot()
moc.write('xmm_MOC_uniq.fits')
table = Table.from_pandas(df)
table.write('xmm_MOC_radec.fits',format='fits',overwrite=True)
%ls -lh
Explanation: MOC catalogs now for LaMassa's XMM
End of explanation
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
res = radec_2_moc('Chandra_multiwavelength.fits','RA','DEC','RADEC_ERR')
print res
res = radec_2_moc('XMM_multiwavelength_cat.fits','RA','DEC','RADEC_ERR')
print res
%ls -lh
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()
print_fits_columns('photometry/hers/hers_catalogue_3sig250_no_extended.fits')
res = radec_2_moc('photometry/hers/hers_catalogue_3sig250_no_extended.fits','RA','DEC')
print res
print_fits_columns('photometry/galex/S82_gmsc_chbrandt.fit')
res = radec_2_moc('photometry/galex/S82_gmsc_chbrandt.fit','ra','dec','poserr')
print res
print_fits_columns('photometry/sdss/Stripe82_photo_chbrandt.fit')
res = radec_2_moc('photometry/sdss/Stripe82_photo_chbrandt.fit','ra','dec')
print res
print_fits_columns('photometry/shela/shela_stripe82_v1.3_cat.fits')
res = radec_2_moc('photometry/shela/shela_stripe82_v1.3_cat.fits','SDSS_RA','SDSS_DEC')
print res
print_fits_columns('photometry/spies/SpIES_ch1ch2_allaor_5s_bothchan_final.fits')
res = radec_2_moc('photometry/spies/SpIES_ch1ch2_allaor_5s_bothchan_final.fits','RA','DEC')
print res
print_fits_columns('photometry/unwise/brandt.fits')
res = radec_2_moc('photometry/unwise/brandt.fits','ra','dec')
print res
print_fits_columns('photometry/vla/first_14dec17.fits')
res = radec_2_moc('photometry/vla/first_14dec17.fits','RA','DEC')
print res
Explanation: Automating MOC generation from fits files
We need, to generate a coverage map, the following parameters:
* R.A. vector
* Dec. vector
* size of the elements
Object::
    given params:
        _ra     -&gt; &lt;vector&gt;
        _dec    -&gt; &lt;vector&gt;
        _radius -&gt; &lt;vector&gt;
    computed params:
        _healpix-level
        _healpix-nside
    output:
        _moc-uniq &lt;- &lt;filename&gt;
        _moc-obj  &lt;- &lt;filename&gt;
End of explanation
<END_TASK>