<jupyter_start><jupyter_text>[healpix]: https://healpix.jpl.nasa.gov[swift]: http://swift.asdc.asi.it/[swift_deepsky]: https://github.com/chbrandt/swift_deepsky Sky discretization with HEALpixThe task we want to accomplish is that of covering a given region of the sky with circles. Not a new problem to the human existence, I want to do this exercise using [Healpix][healpix].Healpix will be used here for binning a series of coordinates; the coordinates may be a real list of sources or an artifical one.In any case, the binning will downsample the list of coordinates.An application for the downsampling is the definition of observational pointings on a survey, where we have to define the minimum amount of to cover a specific region.The pointings are defined by a central coordinate -- RA,Dec -- and a radius value defining the size of the circle (i.e, field of view) each pointing covers. Swift/XRTAt the end of the day we want to define the coverage of a region as if done by the [Swift/XRT][swift] telescope.Such list of pointings will then be used by our [Swift-DeepSky][swift_deepsky] pipeline.The Swift' field-of-view is of $12\ arcmin$.<jupyter_code>def naive_region2pointings(ramin, ramax, decmin, decmax, radius):
    from astropy import units

    bbox = dict(ramin=ramin * units.degree,
                ramax=ramax * units.degree,
                decmin=decmin * units.degree,
                decmax=decmax * units.degree)

    try:
        rad = radius.to('arcmin')
    except:
        rad = radius * units.arcmin

    # get the closest smaller size of healpix elements
    from moc import utils
    level = utils.size_to_level(rad, truncate=True)
    #
    from moc import core
    dsize = core.HEALPIX_LEVELS[level]

    # create a grid of (fake) coordinates, to then create a MOC table from it
    step_size = dsize.to('deg').value/2**0.5

    import numpy as np
    ra_vec = np.arange(bbox['ramin'].value, bbox['ramax'].value, step_size)
    dec_vec = np.arange(bbox['decmin'].value, bbox['decmax'].value, step_size)
    
    import itertools
    grid = np.asarray(list(itertools.product(ra_vec,dec_vec))).T
    
    return grid
from astropy import units
radius = 12 * units.arcmin
%time grid = naive_region2pointings(25, 30, 10, 15, radius)

ra_vec,dec_vec = grid
print(len(ra_vec))
# from astropy.coordinates import SkyCoord
# coords = SkyCoord(grid,unit='deg')<jupyter_output><empty_output><jupyter_text>A list of fixed-order Healpix elements<jupyter_code>def coords2hpix(ra, dec, level):
    """
    Input:
     - radec : array shape = (N,2)
     - level : healpix level representing the size of each element
    """
    from healpy.pixelfunc import ang2pix
    
    nside = 2**level

#     pix_to_visit = []
#     for i in range(len(ra)):
#         ipix = ang2pix(nside, ra[i], dec[i], nest=True, lonlat=True)
#         pix_to_visit.append(ipix)

    def a2p(radec):
        return ang2pix(nside, radec[0], radec[1], nest=True, lonlat=True)
        
#     pix_to_visit = list(map(a2p, zip(ra,dec)))
    pix_to_visit = map(a2p, zip(ra,dec))
        
    return pix_to_visit


def hpix2coords(hpixs, level):
    """
    Transform healpix elements to coordinates
    
    Input:
     - hpixs : list of healpix elements
     - level : healpix level to consider
    """
    from healpy.pixelfunc import pix2ang

    nside = 2**level

#     coords_to_visit = []
#     for ipix in hpixs:
#         c = pixelfunc.pix2ang(nside, ipix, nest=True, lonlat=True)
#         coords_to_visit.append(c)
    def p2a(ipix):
        return pix2ang(nside, ipix, nest=True, lonlat=True)
    
#     coords_to_visit = list(map(p2a, hpixs))
    coords_to_visit = map(p2a, hpixs)
    
    return coords_to_visit


def coords_binning_hpix(ra, dec, level, unique=False):
    """
    Return the healpix binning of 'ra,dec' at some 'level'
    
    Input:
     - radec : array shape = (N,2)
     - level : healpix level representing the size of each element
     - unique: if True, remove (binned) duplicates
    """
    
    pix_to_visit = coords2hpix(ra, dec, level)
    
    if unique:
        pix_to_visit = set(pix_to_visit)
        
    return hpix2coords(pix_to_visit, level)
    

def coords_binning(ra, dec, rad, unique=False, truncate=False):
    """
    Return list of binned ra,dec in closest healpix size to 'rad'
    
    Input:
     - ra       : list of (degree) RA coordinates
     - dec      : list of (degree) Dec coordinates
     - rad      : radius/size value to bin;
                  the actual value will be defined by one of Healpix levels;
                  the closest level may be above or below, see 'truncate'
     - unique   : if True, remove duplicated binned coordinates;
                  if False, returned list matches input ra,dec
     - truncate : if True, use the larger-closest healpix element size;
                  if False, use the smaller-closest healpix element
    """
    assert len(ra) == len(dec)
    assert rad > 0
    
    from moc import utils
    level = utils.size_to_level(rad, truncate=truncate)
    
    coords_to_visit = coords_binning_hpix(ra, dec, level, unique=unique)
    
    return coords_to_visit,level

    ## We don't need the neighbours, but if that was the case:
    #
    #nipixs = pixelfunc.get_all_neighbours(nside, c[0], c[1], nest=True, lonlat=True)
    #cns = []
    #for nipix in nipixs:
    #    cns.append(pixelfunc.pix2ang(nside, nipix, nest=True, lonlat=True))
    #coords_to_visit.extend(cns)
        
%time coords_to_visit,level = coords_binning(ra_vec, dec_vec, radius, True, False)
from moc import core
x_hp,y_hp = list(zip(*coords_to_visit))
r_hp = core.HEALPIX_LEVELS[level].to('deg').value
r_og = radius.to('deg').value

r_middle = (r_hp + r_og)/2

print('radii:', r_hp, r_og, r_middle)
from bokeh.plotting import figure
from bokeh.io import output_notebook, show

output_notebook()

fig = figure()
fig.circle(x_hp,y_hp,radius=r_hp, fill_alpha=0.1, color='red')
# fig.circle(x_hp,y_hp,radius=r_hp, alpha=0.1)
# fig.circle(x_hp,y_hp,radius=r_middle, fill_alpha=0, line_color='black')
show(fig)<jupyter_output><empty_output><jupyter_text>A real catalogLet's build the pointings visualization from a real list of RA,DecWe'll take the LaMassa catalog, in the file `lamassa_chandracsv`<jupyter_code>import pandas

lamassa = pandas.read_csv('lamassa_chandra.csv')
lamassa.describe()
from astropy import units
ra_vec = lamassa['RAdeg']
dec_vec = lamassa['DEdeg']
radius = 12 * units.arcmin

coords_to_visit,level = coords_binning(ra_vec, dec_vec, radius, True, False)
x,y = list(zip(*coords_to_visit))
r = radius.to('deg').value
print(len(x))
from bokeh.plotting import figure
from bokeh.io import output_notebook, show

output_notebook()

fig = figure(x_range=(0,30))
fig.circle(x, y, radius=r, fill_alpha=1, color='blue')
fig.circle(ra_vec, dec_vec, radius=r, fill_alpha=0.1, color='red')
show(fig)<jupyter_output><empty_output><jupyter_text>Create the list of pointings<jupyter_code>import pandas

df = pandas.DataFrame({'ra':x, 'dec':y, 'radius':radius.to('arcmin').value})
df
# df = df.sample(frac=1).reset_index(drop=True)
def sds_cmdline(df_row):
    cmdline = 'docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky'
    ra = df_row['ra']
    dec = df_row['dec']
    radius = df_row['radius']
    lbl = '{:d}_{:d}_{:d}'.format(int(df_row['index']), int(ra), int(dec))
    cmdline += ' --ra {:.5f} --dec {:.5f} --radius {:.3f} --label {!s}'
    cmdline = cmdline.format(ra, dec, radius, lbl)
#     print(cmdline)
    return cmdline
    
cmdlines = df.reset_index().apply(sds_cmdline, axis=1)
cmdlines.to_csv('swift_deepsky_pointings_list.txt', header=False, index=False)
%cat swift_deepsky_pointings_list.txt<jupyter_output>docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 0.43945 --dec 0.74606 --radius 12.000 --label 0_0_0
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 359.38477 --dec 0.59684 --radius 12.000 --label 1_359_0
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 359.47266 --dec 0.67145 --radius 12.000 --label 2_359_0
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 349.54102 --dec 0.74606 --radius 12.000 --label 3_349_0
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 349.45312 --dec 0.82067 --radius 12.000 --label 4_349_0
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 0.61523 --dec 0.89528 --radius 12.000 --label 5_0_0
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 0.61523 --dec 1.04451 --radius 12.000 --label 6_0_1
docker run --rm -v /scratch/work:/work chbrandt/swift_deepsky --ra 359.29688 --dec 0.67145 --radius 12.000 --label 7_359_0
docker run --rm -v /[...]
{'fraction_non_alphanumeric': 0.11824750544787246, 'fraction_numerical': 0.0277554765454754, 'mean_word_length': 3.4354018311291963, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 17, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '276510', 'license': 'MIT', 'n_tokens_mistral': 3308, 'n_tokens_neox': 3081, 'n_words': 838, 'path': 'region_pointings_healpix_realcatalog.ipynb', 'repo_name': 'chbrandt/moc', 'size': '8719'}