MOCPy introduction

MOCPy is a python library for creating, manipulate and parse MOCs (Multi-Order Coverage maps). A MOC describes any arbitrary region on the sky. MOCs can be used to:

  • Represent the spatial footprint of a catalog (source and/or image survey).

  • Compare the footprints, perform fast intersections, unions, differences.

  • Filter an astropy table by discarding all the sources that do not lie in the MOC region.

MOCPy’s code can be found on GitHub: https://github.com/cds-astro/mocpy You can install it: pip install mocpy

MOCPy’s documentation: https://cds-astro.github.io/mocpy/

[1]:

import astropy.units as u
import matplotlib.pyplot as plt
import mocpy
from astropy.coordinates import Angle, SkyCoord
from astropy.visualization.wcsaxes.frame import EllipticalFrame
from astroquery.cds import cds
from astroquery.simbad import Simbad
from astroquery.vizier import Vizier
from ipyaladin import Aladin
from ipywidgets import Box, Layout, widgets
from matplotlib import patheffects
from mocpy import WCS
from regions import CircleSkyRegion

mocpy.__version__

[1]:

'0.13.0'

## Use astroquery.cds to get spatial footprints (MOCs)

MOCs can be retrieved from astroquery.cds. astroquery.cds offers a Python access API to the MOCServer that stores ~20000 metadata and MOCs of Vizier catalogues and ~500 metadata and MOCs of HiPS image surveys.

astroquery.cds documentation https://astroquery.readthedocs.io/en/latest/cds/cds.html#getting-started

Let’s retrieve:

  • The MOC representing the footprint of all the HST combined surveys (see the astroquery.cds documentation, an example is given about that) at the order 8 (i.e. the precision of the MOC, ~13 arcmin)

  • The MOC representing the footprint of SDSS9: ID=’CDS/P/SDSS9/color’

[2]:

# HST MOC footprint
# -----------------
# We want to retrieve all the HST surveys i.e. the HST surveys covering any region of the sky.
allsky = CircleSkyRegion(SkyCoord(0, 0, unit="deg"), Angle(180, unit="deg"))
hst_moc = cds.query_region(
    region=allsky,
    # We want a MOCpy object instead of an astropy table
    return_moc=True,
    # The order of the MOC
    max_norder=11,
    # Expression on the ID meta-data
    meta_data="ID=*HST*",
)

# SDSS9
# -----
sdss_moc = cds.find_datasets(meta_data="ID=CDS/P/SDSS9/color", return_moc=True)

# GALEX
# -----
# galex = cds.find_datasets(meta_data="ID=CDS/P/GALEXGR6/AIS/color", return_moc=True)

[3]:

type(sdss_moc)

[3]:

mocpy.moc.moc.MOC

## Manipulate MOCs using MOCPy

astroquery.cds returns mocpy.MOC typed objects. Use MOCPy (see the API of the mocpy.MOC class https://cds-astro.github.io/mocpy/stubs/mocpy.MOC.html#mocpy.MOC) to manipulate them, for example you could:

  • Compute their intersection/union

  • Serialize them to FITS/json, save them to FITS files

  • Filter an astropy table to keep only the sources being on a MOC (the intersection between sdss and the hst surveys).

[4]:

%%time
sdss_and_hst_moc = sdss_moc.intersection(hst_moc)
sdss_moc.serialize(format="fits")
sdss_moc.save("sdss_moc.fits", format="fits", overwrite=True)

CPU times: user 6.07 ms, sys: 8.07 ms, total: 14.1 ms
Wall time: 13.4 ms

## Plot a MOC using matplotlib

Let’s see how to plot a MOC using matplotlib. There is an example of that on the MOCPy’s documentation: https://cds-astro.github.io/mocpy/examples/examples.html#loading-and-plotting-the-moc-of-sdss

We use matplotlib andMOCPy to draw the MOCs of HST and SDSS that we downloaded from astroquery.cds.

MOCPy offers an interface to create a WCS:

MOCPy offers 2 methods taking a matplotlib.axe.Axe and drawing into it either:

  • the full area of the MOC (mocpy.MOC.fill)

  • only its perimeter (mocpy.MOC.border)

These methods accept additional stylistic kwargs defined by matplotlib: https://matplotlib.org/api/_as_gen/matplotlib.patches.PathPatch.html#matplotlib.patches.PathPatch

[5]:

fig = plt.figure(figsize=(20, 10))

with WCS(
    fig,
    fov=160 * u.deg,
    center=SkyCoord(0, 0, unit="deg", frame="icrs"),
    coordsys="icrs",
    rotation=Angle(0, u.degree),
    projection="AIT",
) as wcs:
    # Create a matplotlib axe and give it a astropy.wcs.WCS-like object
    ax = fig.add_subplot(1, 1, 1, projection=wcs, frame_class=EllipticalFrame)

    # Fill the SDSS MOC in red with an opacity of 70%
    sdss_moc.fill(
        ax=ax,
        wcs=wcs,
        edgecolor="r",
        facecolor="r",
        linewidth=0,
        fill=True,
        alpha=0.7,
        label="sdss9 footprint",
    )
    # Draw its perimeter in black
    sdss_moc.border(ax=ax, wcs=wcs, color="black", alpha=0.5)

    # Fill the HST surveys MOC in green with an opacity of 70%
    hst_moc.fill(
        ax=ax,
        wcs=wcs,
        edgecolor="g",
        facecolor="g",
        linewidth=0,
        fill=True,
        alpha=0.7,
        label="hst all surveys footprint",
    )
    # Draw its perimeter in black
    hst_moc.border(ax=ax, wcs=wcs, color="black", alpha=0.5)

    # Fill the intersection MOC in blue
    sdss_and_hst_moc.fill(
        ax=ax,
        wcs=wcs,
        edgecolor="b",
        facecolor="b",
        linewidth=0,
        fill=True,
        alpha=0.7,
        label="intersection",
    )
    # Draw its perimeter in black
    sdss_and_hst_moc.border(ax=ax, wcs=wcs, color="black", alpha=0.5)

# Usual matplotlib calls
plt.title("Using matplotlib to vizualize MOCs")
plt.xlabel("ra")
plt.ylabel("dec")
plt.legend()
plt.grid(color="black", linestyle="dotted")
path_effects = [patheffects.withStroke(linewidth=3, foreground="black")]
ax.coords[0].set_ticklabel(color="white", path_effects=path_effects)
plt.show()
plt.close()
../../_images/_collections_notebooks_00-MOCpy_introduction_11_0.png

## Filter an astropy.Table by a MOC

  1. Retrieve a catalog table from Vizier (e.g. II/50). Add the columns ‘_RAJ2000’ and ‘_DEJ2000’ to the outputs. MOCPy needs the positions for filtering the table.

  2. Filter the table to get only the sources that lie into intersection MOC.

[ ]:

viz = Vizier(columns=["*", "_RAJ2000", "_DEJ2000"])
viz.ROW_LIMIT = -1
# Photometric standard stars (tables II and IV of paper)
tables = viz.get_catalogs("II/50")
tables

[ ]:

table = tables[0]
table

[ ]:

idx_inside = sdss_and_hst_moc.contains_lonlat(
    table["_RAJ2000"].T * u.deg,
    table["_DEJ2000"].T * u.deg,
)
sources_inside = table[idx_inside]
sources_inside

## Run Aladin-Lite inside a jupyter notebook: ipyaladin

Aladin-Lite can be embedded into a jupyter notebook: Follow the readme on GitHub for installing it: https://github.com/cds-astro/ipyaladin

[ ]:

aladin = Aladin()
aladin

[ ]:

aladin.target = "messier 51"
aladin.fov = 1
aladin.height = 600

[ ]:

aladin.coo_frame = "icrs"

Change the image survey, go to https://aladin.unistra.fr/hips/list (Part 1. HiPS sky maps) to test with different image surveys! A list of good HiPS I like for testing: - P/2MASS/color - P/PanSTARRS/DR1/color-z-zg-g - P/SPITZER/color - P/SDSS9/color - P/GALEXGR6/AIS/color - P/Mellinger/color

[ ]:

aladin.survey = "P/SDSS9/color"

Add a MOC in the aladin-lite view!

[ ]:

# First use MOCPy to serialize our MOC into the JSON format
moc_serialized_json = sdss_and_hst_moc.serialize("json")
aladin.add_moc_from_dict(
    moc_serialized_json,
    {"color": "blue", "opacity": 0.3, "fill": True, "edge": False, "perimeter": True},
)

[ ]:

# Add astropy source tables to the aladin lite viewer
aladin.add_table(sources_inside)
aladin.add_table(table[~idx_inside])

[ ]:

# change the fov and target
aladin.target = "13 04 4.193 -03 34 13.54"
aladin.fov = 11

Pass Python callback functions to Javascript event listeners

[ ]:

aladin = Aladin(target="M 1", fov=0.2)
aladin

[ ]:

table = Simbad.query_region("M 1", radius=0.1 * u.deg)
aladin.add_table(table)
aladin.height = 600

[ ]:

def get_ra_dec_data(source):
    """Return the MainID, and coordinates of the clicked sources."""
    return [source["data"]["MAIN_ID"], source["ra"], source["dec"]]


def hover_source_callback(source):
    """Return data corresponding to the hovered sources."""
    return source["data"]


# When trigerred, the listeners on the js side of the application will send a
# json object whose parameter data
# will be used by the python functions
# (data is a litteral object on the js side, it will be converted as a
# dictionary object on the Python side)
aladin.add_listener("objectHovered", hover_source_callback)
aladin.add_listener("objectClicked", get_ra_dec_data)

Advanced example

Create a rectangular selection using ipywidgets and gives back the selected sources to the user

[ ]:

aladin = Aladin(layout=Layout(width="70%"), target="M 1", fov=0.2)

button = widgets.Button(description="Select")


def on_button_clicked(b):  # noqa: ARG001
    """Trigger the rectangular selection tool."""
    aladin.rectangular_selection()


button.on_click(on_button_clicked)
table_info = widgets.HTML(layout=Layout(height="auto", overflow="auto"))

box_layout = Layout(
    display="flex",
    flex_flow="row",
    align_items="stretch",
    width="100%",
    overflow="hidden",
    height="100vh",
    margin="-100px 0 0 0",
    padding="100px 0 0 0 ",
)
box = Box(children=[aladin, button, table_info], layout=box_layout)
box

[ ]:

table = Simbad.query_region("M 1", radius=0.1 * u.deg)
aladin.add_table(table)


def process_result(sources):
    """Output an html table."""
    s = '<table border="1">'
    s += "<tr><th>MAIN_ID</th><th>RA</th><th>DEC</th></tr>"
    for source in sources:
        s += "<tr><td>{}</td><td>{}</td><td>{}</td></tr>".format(
            source["data"]["MAIN_ID"],
            source["data"]["RA"],
            source["data"]["DEC"],
        )
    s += "</table>"
    table_info.value = s


aladin.add_listener("select", process_result)