In [1]:
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[1]:
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:
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:
In [2]:
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__
Download ReadMe and chandra.dat files and save them inside ./data/ dir
In [3]:
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)
In [4]:
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/'
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.
In [5]:
from astropy.table import Table
chandra = Table.read('data/chandra.dat',readme='data/ReadMe',format='ascii.cds')
In [6]:
chandra # Notice the '-999' values
Out[6]:
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:
In [7]:
from astropy.table import Table
chandra = Table.read('data/chandra.dat',readme='data/ReadMe_fix',format='ascii.cds')
chandra
Out[7]:
That's it. After declaring the null values with all the significant numbers following the Format, such values are properly handled.
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.
In [8]:
# 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)
In [9]:
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)
In [10]:
radec = zip(chandra['RAdeg'],chandra['DEdeg'])
hpix = [ healpix_radec2pix(nside,ra,dec) for ra,dec in radec ]
Here it is, the MOC catalog (the list of elements to be more precise):
In [11]:
hpix
Out[11]:
The plot made by MOCPy:
In [12]:
moc = mocpy.MOC()
moc.add_pix_list(level,hpix)
moc.plot()
moc.write('data/MOC_chandra.fits')
And here after importing it to Aladin:
In [13]:
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>
''')
Out[13]:
In [ ]: