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 astropy.wcs import WCS
from astroquery.mocserver import MOCServer
from astroquery.vizier import Vizier
from ipyaladin import Aladin
from matplotlib import patheffects
from regions import CircleSkyRegion

mocpy.__version__

[1]:

'0.15.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 = MOCServer.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 = MOCServer.find_datasets(
    meta_data="ID=CDS/P/SDSS9/color",
    return_moc=True,
    max_norder=11,
)

[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 12.1 ms, sys: 0 ns, total: 12.1 ms
Wall time: 11.6 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))

wcs = WCS(
    {
        "naxis": 2,
        "naxis1": 1620,
        "naxis2": 810,
        "crpix1": 810.5,
        "crpix2": 405.5,
        "cdelt1": -0.2,
        "cdelt2": 0.2,
        "ctype1": "RA---AIT",
        "ctype2": "DEC--AIT",
    },
)

# Create a matplotlib axe and give it a astropy.wcs.WCS 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)

# we set the axis off by half a pixel otherwise they are slightly outside the sphere
# see astropy issue for the half-pixel trick https://github.com/astropy/astropy/issues/10201
ax.set_xlim(-0.5, 1620 - 0.5)
ax.set_ylim(-0.5, 810 - 0.5)
ax.set_aspect(1.0)

# 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_II50 = tables[0]
table_II50

[ ]:

idx_inside = sdss_and_hst_moc.contains_lonlat(
    table_II50["_RAJ2000"].T * u.deg,
    table_II50["_DEJ2000"].T * u.deg,
)
sources_inside = table_II50[idx_inside]
sources_outside = table_II50[~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!

[ ]:

aladin.add_moc(
    sdss_and_hst_moc,
    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, color="green", shape="rhomb")
aladin.add_table(sources_outside, color="red", shape="cross")

[ ]:

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