<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Building a MOC from a CDS/Vizier catalog
After a couple of conversations I had during the week of the 2nd ASTERICS VO school I've got two small TODO's hanging on my list, this is an answer to both. First one (to François) is an observation about CDS/Vizier LaMassa 2016 catalog's metadata; the second one (to Markus) is about MOC catalogs format. I will not get into the details -- since the respective discussion had already being done --, it is sufficient to say that
Step1:  Download ReadMe and chandra.dat files and save them inside ./data/ dir
Step2:  Dealing with null values from Vizier metadata
The goal here, as said before, is to show the "bug" in Vizier ReadMe (description) file when dealing with null values not properly formatted.
We start by opening the chandra table and noticing the values -999 not being properly handled by Astropy as null values.
Step3:  We can see records in columns logLSoft, logLHard, logLFull, for example, showing the values -999.0.
Although we already suspect, we go to ReadMe and double-check it to see what is there about Null values.
For those three example columns we see ?=-999, which is the right --although truncated/integer-- value.
For some reason, Astropy is not handling it as it should.
To fix this, we have sync the number of significat digits of the null values with the (declared) format.
For instance, those columns have a Format=F7.2 and so should the null values be declared as ?=-999.00.
Changing such signatures for those columns (logLSoft, logLHard and logLFull) and saving them to file ReadMe_fix give us the following
Step7:  That's it. After declaring the null values with all the significant numbers following the Format, such values are properly handled.
Generating a MOC catalog
Now we go through a MOC catalog creation. There is no problem here, just to answer Markus what I understand what a MOC is
Step8:  Here it is, the MOC catalog (the list of elements to be more precise)
Step9:  The plot made by MOCPy
Step10:  And here after importing it to Aladin
<ASSISTANT_TASK:>
Python Code:
baseurl = 'ftp://cdsarc.u-strasbg.fr/pub/cats/J/ApJ/817/172/'
readme_file = 'ReadMe'
chandra_file = 'chandra.dat'
import astropy
print "astropy version:",astropy.__version__
import mocpy
print "mocpy version:",mocpy.__version__
import healpy
print "healpy version:",healpy.__version__
Explanation: Building a MOC from a CDS/Vizier catalog
After a couple of conversations I had during the week of the 2nd ASTERICS VO school I've got two small TODO's hanging on my list, this is an answer to both. First one (to François) is an observation about CDS/Vizier LaMassa 2016 catalog's metadata; the second one (to Markus) is about MOC catalogs format. I will not get into the details -- since the respective discussion had already being done --, it is sufficient to say that:
some of Lamassa's null value are not properly read when using Astropy;
what I understand about a MOC catalog.
The data to be used is provided by Vizier, the LaMassa et al, 2016, ApJ, 817, 172, in particular the ReadMe and chandra.dat files.
Software to handle the catalog will be Astropy and Thomas Boch's mocpy and Healpy.
TOC:
 * Dealing with null values from Vizier metadata
 * Generating a MOC catalog
End of explanation
def download(path,filename,outdir):
    from urllib2 import urlopen
    url = path+filename
    f = urlopen(url)
    data = f.read()
    with open(outdir+filename, "wb") as fp:
        fp.write(data)
import os
if not os.path.isdir('data'):
    os.mkdir('data')
download(baseurl,readme_file,outdir='data/')
download(baseurl,chandra_file,outdir='data/')
!ls 'data/'
Explanation: Download ReadMe and chandra.dat files and save them inside ./data/ dir
End of explanation
from astropy.table import Table
chandra = Table.read('data/chandra.dat',readme='data/ReadMe',format='ascii.cds')
chandra # Notice the '-999' values
Explanation: Dealing with null values from Vizier metadata
The goal here, as said before, is to show the "bug" in Vizier ReadMe (description) file when dealing with null values not properly formatted.
We start by opening the chandra table and noticing the values -999 not being properly handled by Astropy as null values.
End of explanation
from astropy.table import Table
chandra = Table.read('data/chandra.dat',readme='data/ReadMe_fix',format='ascii.cds')
chandra
Explanation: We can see records in columns logLSoft, logLHard, logLFull, for example, showing the values -999.0.
Although we already suspect, we go to ReadMe and double-check it to see what is there about Null values.
For those three example columns we see ?=-999, which is the right --although truncated/integer-- value.
For some reason, Astropy is not handling it as it should.
To fix this, we have sync the number of significat digits of the null values with the (declared) format.
For instance, those columns have a Format=F7.2 and so should the null values be declared as ?=-999.00.
Changing such signatures for those columns (logLSoft, logLHard and logLFull) and saving them to file ReadMe_fix give us the following:
End of explanation
# 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'.
    
    # 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
    }
    
    assert size.unit
    ko = None
    for k,v in healpix_levels.iteritems():
        if v < 2 * size:  # extrapolating the error by one order of magnitude
            break
        ko = k
    return ko
import numpy as np
from astropy import units as u
median_positional_error = np.median(chandra['e_Pos']) * u.arcsec
level = size2level(median_positional_error)
nside = 2**level
print "Typical (median) position error: \n{}".format(median_positional_error)
print "\nCorrespondig healpix level: {} \n\t and nsize value: {}".format(level,nside)
def healpix_radec2pix(nside, ra, dec, nest=True):
    
    convert ra,dec to healpix elements
    
    def radec2thetaphi(ra,dec):
        
        convert equatorial ra, dec in degrees
        to polar theta, phi in radians
        
        def ra2phi(ra):
            import math
            return math.radians(ra)
        def dec2theta(dec):
            import math
            return math.pi/2 - math.radians(dec)
        _phi = ra2phi(ra)
        _theta = dec2theta(dec)
        return _theta,_phi
    
    import healpy
    _theta,_phi = radec2thetaphi(ra, dec)
    return healpy.ang2pix(nside, _theta, _phi, nest=nest)
radec = zip(chandra['RAdeg'],chandra['DEdeg'])
hpix = [ healpix_radec2pix(nside,ra,dec) for ra,dec in radec ]
Explanation: That's it. After declaring the null values with all the significant numbers following the Format, such values are properly handled.
Generating a MOC catalog
Now we go through a MOC catalog creation. There is no problem here, just to answer Markus what I understand what a MOC is: a list of (unique) element numbers. The section 2.3.1 NUNIQ packing of the IVOA MOC document version-1 explains how the conversion between the two kinds of elements representation.
The lines below will use the catalog we just have in hands, chandra, to build the MOC.
First, I will compute Healpix level/nside values based on the mean positional error of the catalog and then the MOC elements are computed from (RA,Dec) using HealPy.
MOCPy is used at the end to plot the MOC elements; Also Aladin is used to have a better view of the elements.
Files can be download from data/ directory.
End of explanation
hpix
Explanation: Here it is, the MOC catalog (the list of elements to be more precise):
End of explanation
moc = mocpy.MOC()
moc.add_pix_list(level,hpix)
moc.plot()
moc.write('data/MOC_chandra.fits')
Explanation: The plot made by MOCPy:
End of explanation
from IPython.display import HTML
HTML('''
<figure>
  <img src="data/MOC_on_Aladin.png" alt="MOC printed on Aladin">
  <figcaption>Figure 1: MOC printed on Aladin</figcaption>
</figure>
''')
Explanation: And here after importing it to Aladin:
End of explanation
<END_TASK>