User documentation¶

SMOC (a.k.a MOC): Spatial coverages¶

Gallery of notebooks examples using SMOCs¶

MOCPy introduction

Compute the borders of a MOC

Create a MOC from a filtered astropy Table

Overlay the MOC of Spitzer on a 2MASS cutout

Create a MOC from an astropy Table

Create a MOC from a polygon

Intersection between MOCs

Construct a MOC from an image with a mask

Get the MOC corresponding to a table

Query the sources of a VizieR table that are in a sky region

Here are some code examples manipulating MOC objects.

Examples use cases¶

Loading and plotting the MOC of SDSS¶

This example:

Load the coverage of SDSS from a FITS file.
Plot the MOC by:

Defining a matplotlib figure.
Defining an astropy WCS representing the field of view of the plot.
Call the mocpy.moc.MOC.fill() and mocpy.moc.MOC.border() so that mocpy plot on a matplotlib axis.
Set the axis labels, a title, enable the grid and plot the final figure.

import astropy.units as u
import matplotlib.pyplot as plt
from astropy.coordinates import Angle

from mocpy import MOC

# Load a MOC
filename = "./../../resources/P-SDSS9-r.fits"
moc = MOC.from_fits(filename)
# Plot the MOC using matplotlib
fig = plt.figure(111, figsize=(15, 10))
# Define a astropy WCS easily
wcs = moc.wcs(fig, coordsys="icrs", rotation=Angle(0, u.degree), projection="AIT")
ax = fig.add_subplot(1, 1, 1, projection=wcs)
# Call fill with a matplotlib axe and the `~astropy.wcs.WCS` wcs object.
moc.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color="green")
moc.border(ax=ax, wcs=wcs, alpha=0.5, color="black")
plt.xlabel("ra")
plt.ylabel("dec")
plt.title("Coverage of P-SDSS9-r")
plt.grid(color="black", linestyle="dotted")
plt.show()

(Source code, png, hires.png, pdf)

Intersection between GALEX and SDSS¶

This example:

Load the coverages of SDSS and GALEX from FITS files.
Compute their intersection
Compute their union
Plot the resulting intersection and union on a same matplotlib axis.

import astropy.units as u
import matplotlib.pyplot as plt
from astropy.coordinates import Angle, SkyCoord

from mocpy import MOC, WCS

# Load Galex and SDSS
sdss = MOC.from_fits("./../../resources/P-SDSS9-r.fits")
galex = MOC.from_fits("./../../resources/P-GALEXGR6-AIS-FUV.fits")
# Compute their intersection
inter = sdss & galex
union = sdss + galex
# Plot the MOC using matplotlib
fig = plt.figure(111, figsize=(10, 10))
# Define a astropy WCS easily
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:
    ax = fig.add_subplot(1, 1, 1, projection=wcs)
    # Call fill with a matplotlib axe and the `~astropy.wcs.WCS` wcs object.
    union.fill(
        ax=ax,
        wcs=wcs,
        alpha=0.5,
        fill=True,
        color="red",
        linewidth=0,
        label="Union",
    )
    union.border(ax=ax, wcs=wcs, alpha=1, color="red")

    inter.fill(
        ax=ax,
        wcs=wcs,
        alpha=0.5,
        fill=True,
        color="green",
        linewidth=0,
        label="Intersection",
    )
    inter.border(ax=ax, wcs=wcs, alpha=1, color="green")
    ax.legend()

plt.xlabel("ra")
plt.ylabel("dec")
plt.title("Logical operations between SDSS and GALEX")
plt.grid(color="black", linestyle="dotted")
plt.show()

(Source code, png, hires.png, pdf)

Create a MOC from a concave polygon¶

This example shows how to call mocpy.moc.MOC.from_polygon() or mocpy.moc.MOC.from_polygon_skycoord().

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import Angle, SkyCoord
from mocpy import MOC, WCS

# The set of points delimiting the polygon in deg
vertices = np.array(
    [
        [18.36490956, 5.0],
        [15.7446692, 6.97214743],
        [16.80755056, 10.29852928],
        [13.36215502, 10.14616136],
        [12.05298691, 13.10706197],
        [9.54793022, 10.4556709],
        [8.7677627, 7.80921809],
        [9.71595962, 5.30855011],
        [7.32238541, 6.44905255],
        [0.807265, 6.53399616],
        [1.08855146, 3.51294225],
        [2.07615384, 0.7118289],
        [3.90690158, -1.61886929],
        [9.03727956, 2.80521847],
        [9.22274427, -4.38008174],
        [10.23563378, 4.06950324],
        [10.79931601, 3.77655576],
        [14.16533992, 1.7579858],
        [19.36243491, 1.78587203],
        [15.31732084, 5.0],
    ],
)
skycoord = SkyCoord(vertices, unit="deg", frame="icrs")
# A point that we say it belongs to the inside of the MOC
inside = SkyCoord(ra=10, dec=5, unit="deg", frame="icrs")
moc = MOC.from_polygon_skycoord(skycoord, max_depth=9)
moc.save("polygon_moc.fits", format="fits", overwrite=True)

# Plot the MOC using matplotlib
fig = plt.figure(111, figsize=(10, 10))
# Define a astropy WCS easily
with WCS(
    fig,
    fov=20 * u.deg,
    center=SkyCoord(10, 5, unit="deg", frame="icrs"),
    coordsys="icrs",
    rotation=Angle(0, u.degree),
    # The gnomonic projection transforms great circles into straight lines.
    projection="TAN",
) as wcs:
    ax = fig.add_subplot(1, 1, 1, projection=wcs)
    # Call fill with a matplotlib axe and the `~astropy.wcs.WCS` wcs object.
    moc.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color="red", linewidth=1)
    moc.border(ax=ax, wcs=wcs, alpha=1, color="red")

plt.xlabel("ra")
plt.ylabel("dec")
plt.title("MOC from a polygon")
plt.grid(color="black", linestyle="dotted")
plt.show()

(Source code, png, hires.png, pdf)

Get the border(s) of a MOC¶

This example shows how to call mocpy.moc.MOC.get_boundaries. The borders are returned as a list of SkyCoord each defining one border. In this example:

The sky coordinates defining the border of the MOC are projected to the pixel image system.
Then, a matplotlib path is defined from the projected vertices.
Finally the path is plot on top of the MOC.

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import Angle, SkyCoord
from matplotlib.patches import PathPatch
from matplotlib.path import Path

from mocpy import MOC, WCS

moc = MOC.from_fits("polygon_moc.fits")
skycoords = moc.get_boundaries()

# Plot the MOC using matplotlib
fig = plt.figure(111, figsize=(10, 10))
# Define a astropy WCS easily
with WCS(
    fig,
    fov=20 * u.deg,
    center=SkyCoord(10, 5, unit="deg", frame="icrs"),
    coordsys="icrs",
    rotation=Angle(0, u.degree),
    # The gnomonic projection transforms great circles into straight lines.
    projection="TAN",
) as wcs:
    ax = fig.add_subplot(1, 1, 1, projection=wcs)
    # Call fill with a matplotlib axe and the `~astropy.wcs.WCS` wcs object.
    moc.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color="red", linewidth=1)
    moc.border(ax=ax, wcs=wcs, alpha=1, color="red")

    # Plot the border
    from astropy.wcs.utils import skycoord_to_pixel

    x, y = skycoord_to_pixel(skycoords[0], wcs)
    p = Path(np.vstack((x, y)).T)
    patch = PathPatch(p, color="black", fill=False, alpha=0.75, lw=2)
    ax.add_patch(patch)

plt.xlabel("ra")
plt.ylabel("dec")
plt.title("Get the border(s) of a MOC")
plt.grid(color="black", linestyle="dotted")
plt.show()

(Source code, png, hires.png, pdf)

Gravitational Waves MOCs¶

This example shows the probability confidence regions of gravitational waves. HEALPix cells are given under the uniq pixel notation. Each pixel is associated with a specific value. We can create a MOC from which a GW has x% of change of being localized in it. By definition the MOC which has 100% of chance of containing a GW is the full sky MOC.

import astropy.units as u
import astropy_healpix as ah
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import Angle, SkyCoord
from astropy.io import fits

from mocpy import MOC, WCS

fits_image_filename = "./../../resources/bayestar.multiorder.fits"

max_order = None
with fits.open(fits_image_filename) as hdul:
    hdul.info()
    data = hdul[1].data
    max_order = hdul[1].header["MOCORDER"]

uniq = data["UNIQ"]
probdensity = data["PROBDENSITY"]


level, ipix = ah.uniq_to_level_ipix(uniq)
area = ah.nside_to_pixel_area(ah.level_to_nside(level)).to_value(u.steradian)

prob = probdensity * area


cumul_to = np.linspace(0.5, 0.9, 5)[::-1]
colors = ["blue", "green", "yellow", "orange", "red"]
mocs = [
    MOC.from_valued_healpix_cells(uniq, prob, max_order, cumul_to=c) for c in cumul_to
]


# Plot the MOC using matplotlib


fig = plt.figure(111, figsize=(15, 10))
# Define a astropy WCS easily
with WCS(
    fig,
    fov=50 * u.deg,
    center=SkyCoord(315, 15, unit="deg", frame="icrs"),
    coordsys="icrs",
    rotation=Angle(0, u.degree),
    projection="AIT",
) as wcs:
    ax = fig.add_subplot(1, 1, 1, projection=wcs)
    # Call fill with a matplotlib axe and the `~astropy.wcs.WCS` wcs object.
    for moc, c, col in zip(mocs, cumul_to, colors):
        moc.fill(
            ax=ax,
            wcs=wcs,
            alpha=0.5,
            linewidth=0,
            fill=True,
            color=col,
            label="confidence probability " + str(round(c * 100)) + "%",
        )
        moc.border(ax=ax, wcs=wcs, alpha=0.5, color=col)

    ax.legend()

plt.xlabel("ra")
plt.ylabel("dec")
plt.title("Bayestar")
plt.grid(color="black", linestyle="dotted")
plt.show()

(Source code, png, hires.png, pdf)

Performing computation on the pixels of an FITS image lying in a MOC¶

This example shows how a MOC can filter pixels from a specific FITS image (i.e. associated with a WCS). These pixels can then be retrieved from the image for performing some computations on them: e.g. mean, variance analysis thanks to numpy/scikit-learn…

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.visualization import simple_norm
from astropy.wcs import WCS as WCS

from mocpy import MOC

# load 2MASS cutout covering the galactic plane
hdu = fits.open(
    "http://alasky.u-strasbg.fr/hips-image-services/hips2fits?hips=CDS%2FP%2F2MASS%2FK&width=1200&height=700&fov=30&projection=TAN&coordsys=galactic&rotation_angle=0.0&object=gal%20center&format=fits",
)

# load Spitzer MOC
moc = MOC.from_fits("http://skies.esac.esa.int/Spitzer/IRAC1_bright_ISM/Moc.fits")

# create WCS from 2MASS image header
twomass_wcs = WCS(header=hdu[0].header)

# compute skycoords for every pixel of the image
width = hdu[0].header["NAXIS1"]
height = hdu[0].header["NAXIS2"]

xv, yv = np.meshgrid(np.arange(0, width), np.arange(0, height))
skycoords = twomass_wcs.pixel_to_world(xv, yv)

ra, dec = skycoords.icrs.ra.deg, skycoords.icrs.dec.deg
# Get the skycoord lying in the MOC mask
mask_in_moc = moc.contains(ra * u.deg, dec * u.deg)

img = hdu[0].data
img_test = img.copy()

# Set to zero the pixels that do not belong to the MOC
img_test[~mask_in_moc] = 0

# Plot the filtered pixels image
fig = plt.figure(111, figsize=(15, 10))
ax = fig.add_subplot(1, 1, 1, projection=twomass_wcs)

im = ax.imshow(
    img_test,
    origin="lower",
    norm=simple_norm(hdu[0].data, "sqrt", min_cut=-1, max_cut=150),
)
plt.show()

(Source code, png, hires.png, pdf)

../_images/filter_image_pixels_00_00.png

# Sum of pixels lying in the MOC
np.sum(img[mask_in_moc])

# Sum of pixels of the whole image
np.sum(img)

TMOC: Temporal coverages¶

The TimeMOC class represents a temporal coverage.

Gallery of notebooks examples using TMOCs¶

Loading a TimeMoc from a fits file

STMOC: Space & Time coverages¶

Space-Time coverages are a new feature of mocpy since its version 0.7.0 and are bind spatial and temporal coverages together. The standard description is published by the IVOA here.

Space-Time coverages allow to:

Retrieve the spatial coverage observed by a mission within a set of time frames (i.e. astropy.time.Time ranges).
Retrieve the temporal coverage observed by a mission within a spatial coverage.

As we do for spatial or temporal coverages, one can also perform the union, intersection or difference between two Space-Time coverages.