Merge a1e9d3d into 57cd016

.circleci/config.yml

@@ -8,7 +8,7 @@ jobs:
       - checkout
       - run:
           name: Install dependencies
-          command: /opt/python/cp36-cp36m/bin/pip install numpy scipy astropy pytest-astropy pytest-xdist Cython jinja2 astropy-healpix
+          command: /opt/python/cp36-cp36m/bin/pip install numpy scipy astropy astropy-healpix pytest-astropy pytest-xdist Cython jinja2 spherical-geometry networkx
       - run:
           name: Run tests
           command: PYTHONHASHSEED=42 /opt/python/cp36-cp36m/bin/python setup.py test --parallel=4

.gitignore

@@ -28,6 +28,9 @@ docs/_build
 # Pycharm editor project files
 .idea
 
+# VSCode editor project files
+.vscode/
+
 # Packages/installer info
 *.egg
 *.egg-info

.rtd-environment.yml

@@ -9,6 +9,7 @@ dependencies:
   - numpy
   - cython
   - matplotlib
+  - spherical-geometry
   - sphinx
   - astropy
   - shapely

.travis.yml

@@ -38,7 +38,7 @@ env:
         # Add the current dev version of astropy-healpix waiting for the
         # new release to support Python 3.7.
         # Cython is necessary for pip installing the dev version of astropy-healpix.
-        - PIP_DEPENDENCIES='Cython git+http://github.com/astropy/astropy-healpix.git#egg=astropy-healpix'
+        - PIP_DEPENDENCIES='Cython astropy-healpix networkx spherical-geometry'
         - CONDA_DEPENDENCIES='Cython shapely pytest-arraydiff'
         - CONDA_CHANNELS='astropy'
         - SETUP_XVFB=True

CHANGES.rst

@@ -1,7 +1,11 @@
 0.4 (Unreleased)
 ================
 
-- No changes yet.
+- Add a ``regions.MOCSkyRegion`` and ``regions.MOCPixelRegion`` regions.
+  ``regions.MOCSkyRegion`` allows to define a coverage map of the sky.
+  It relies on what has been done in `mocpy <https://github.com/cds-astro/mocpy/>`__. It
+  offers plotting, logical operations between MOCs features. MOCs can be stored
+  and loaded from FITS files and can be serialized in JSON. [#219]
 
 0.3 (2018-09-09)
 ================

appveyor.yml

@@ -13,7 +13,7 @@ environment:
       # For this package-template, we include examples of Cython modules,
       # so Cython is required for testing. If your package does not include
       # Cython code, you can set PIP_DEPENDENCIES='' and CONDA_DEPENDENCIES=''
-      PIP_DEPENDENCIES: 'Cython git+http://github.com/astropy/astropy-healpix.git#egg=astropy-healpix'
+      PIP_DEPENDENCIES: 'Cython astropy-healpix networkx spherical-geometry'
       CONDA_CHANNELS: "defaults conda-forge astropy"
       CONDA_DEPENDENCIES: "Cython pytest-arraydiff pytest=3.6.4"
       ASTROPY_USE_SYSTEM_PYTEST: 1

docs/index.rst

@@ -32,6 +32,7 @@ User Documentation
    installation
    getting_started
    shapes
+   moc/moc
    contains
    compound
    masks

docs/moc/moc.rst

@@ -0,0 +1,184 @@
+.. _moc:
+
+Multi-Order Coverage maps
+=========================
+
+.. _moc-intro:
+
+Introduction
+------------
+
+A MOC describes the coverage of an arbitrary region on the unit sphere.
+MOCs are usually used for describing the global coverage of catalog/image surveys such as GALEX or SDSS.
+A MOC corresponds to a list of `HEALPix <https://healpix.sourceforge.io/>`__ cells at different depths.
+This class gives you the possibility to:
+
+1. Define `~regions.MOCSkyRegion` objects:
+
+- From a FITS file that stores HEALPix cells (see `~regions.MOCSkyRegion.from_fits`). :ref:`[1] <from_fits>`
+- Directly from a list of HEALPix cells expressed either as a numpy structural array (see `~regions.MOCSkyRegion.from_cells`) or a simple
+    python dictionnary (see `~regions.MOCSkyRegion.from_json`).
+- From a list of sky coordinates (see `~regions.MOCSkyRegion.from_skycoord`, `~regions.MOCSkyRegion.from_lonlat`).
+- From a convex/concave polygon (see `~regions.MOCSkyRegion.from_polygon`). :ref:`[3] <from_polygon>`
+- From a cone (will be implemented in a next version).
+
+2. Perform fast logical operations between `~regions.MOCSkyRegion` objects: :ref:`[4] <logical_operations>`
+
+- The `~regions.MOCSkyRegion.intersection` 
+- The `~regions.MOCSkyRegion.union`
+- The `~regions.MOCSkyRegion.difference`
+- The `~regions.MOCSkyRegion.complement`
+
+3. :ref:`plot_moc`:
+
+- Draw the MOC with its HEALPix cells (see `~regions.MOCSkyRegion.fill`)
+- Draw the perimeter of a MOC (see `~regions.MOCSkyRegion.border`)
+
+4. Get the sky coordinates defining the border(s) of `~regions.MOCSkyRegion` objects (see `~regions.MOCSkyRegion.get_boundaries`).
+
+5. Serialize `~regions.MOCSkyRegion` objects to `astropy.io.fits.HDUList` or JSON dictionary and save it to a file. :ref:`[6] <serialize>`
+
+More informations about MOCs can be found by reading the following `IVOA paper <http://www.ivoa.net/documents/MOC/20140602/REC-MOC-1.0-20140602.pdf>`_.
+
+Examples
+--------
+
+.. _from_fits:
+
+Load a `~regions.MOCSkyRegion` from a FITS file
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+A MOC is often loaded from a FITS file. A FITS file describing a MOC contains an `astropy.io.fits.BinTableHDU`
+table where one can find a list of 64 bits unsigned numbers. Each of these numbers (usually called uniq numbers)
+describes one HEALPix cell at a specific depth.
+
+For the purpose of the example, we will load the MOC describing the coverage
+of the ``P-GALEXGR6-AIS-FUV`` sky survey. It is a rather complex MOC containing lots of holes
+and spread over a good range of different depths.
+
+.. code-block:: python
+
+    >>> from regions import MOCSkyRegion
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> galex = MOCSkyRegion.from_fits(get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions'))
+
+It is also possible to create a MOC from a list of `~astropy.coordinates.SkyCoord`, lon and lat `~astropy.units.Quantity`,
+a python dict storing HEALPix cells etc... Please refer to the `~regions.MOCSkyRegion` API for more details.
+
+
+.. _plot_moc:
+
+Plot a `~regions.MOCSkyRegion`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This example:
+
+- Load the coverage of GALEX from a FITS file.
+- Plot the MOC by:
+
+1. Defining a matplotlib figure.
+2. Defining an astropy WCS representing the field of view of the plot.
+3. Create a `~regions.MOCPixelRegion` from the moc sky region instance and plot it.
+4. Set the axis labels, a title, enable the grid and plot the final figure.
+
+.. plot:: moc/plot_moc.py
+    :include-source: true
+
+It is also possible to use the `regions.MOCSkyRegion.fill` and `regions.MOCSkyRegion.border` methods from `MOCSkyRegion` without casting it to a `MOCPixelRegion`.
+These methods take an `~astropy.wcs.WCS` and a `matplotlib.axes.Axes` to draw the MOC into it.
+
+.. plot:: moc/plot_moc_no_cast.py
+    :include-source: true
+
+.. _from_polygon:
+
+Create `~regions.MOCSkyRegion` from a polygon
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This example creates a `~regions.MOCSkyRegion` instance from a list of sky coordinates.
+These sky coordinates refer to the vertices of a polygon.
+Concave/Convex polygons are accepted but not self intersecting ones.
+Two methods allow you to create a `~regions.MOCSkyRegion` instance from a polygon:
+
+- `~regions.MOCSkyRegion.from_polygon` asks for two `~astropy.units.Quantity` storing the longitudes (resp. the latitudes) of the polygon vertices.
+- `~regions.MOCSkyRegion.from_polygon_skycoord` asks for an `~astropy.coordinates.SkyCoord` storing the vertices of the polygon.
+
+.. plot:: moc/plot_polygon.py
+    :include-source: true
+
+.. _logical_operations:
+
+Intersection between GALEX and SDSS
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This example:
+
+- Load the coverages of SDSS and GALEX from FITS files.
+
+.. code-block:: python
+
+    >>> from regions import MOCSkyRegion
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> galex = MOCSkyRegion.from_fits(get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions'))
+    >>> sdss9 = MOCSkyRegion.from_fits(get_pkg_data_filename('shapes/tests/data/P-SDSS9-r.fits', package='regions'))
+
+- Compute their intersection
+
+.. code-block:: python
+
+    >>> intersection = galex.intersection(sdss9)
+
+- Compute their union
+
+.. code-block:: python
+
+    >>> union = galex.union(sdss9)
+
+- Plot the resulting intersection and union on a same matplotlib axis.
+
+.. plot:: moc/plot_logical_operations.py
+    :include-source: false
+
+Check whether sky coordiantes fall inside a `~regions.MOCSkyRegion`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Given ra and dec `~astropy.units.Quantity` instances, one can check if these sky positions
+lie inside a `~regions.MOCSkyRegion`.
+
+.. code-block:: python
+
+    >>> import numpy as np
+    >>> from astropy import units as u
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> galex = MOCSkyRegion.from_fits(get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions'))
+    >>> # Generate a thousand random positions on the sky 
+    ... ra = np.random.randint(low=0, high=360, size=1000) * u.deg
+    ... dec = np.random.randint(low=-90, high=90, size=1000) * u.deg
+    >>> inside_mask = galex.contains(ra, dec)
+
+.. plot:: moc/plot_contains.py
+    :include-source: false
+
+.. _serialize:
+
+Write a `~regions.MOCSkyRegion` to a FITS/JSON file
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+One can serialize a `~regions.MOCSkyRegion` into a `~astropy.io.fits.hdu.hdulist.HDUList` object or a
+simple python dictionary (JSON format). The following block of code 
+serializes the galex `~regions.MOCSkyRegion` instance defined in :ref:`[1] <from_fits>` to a FITS hdulist.
+
+.. code-block:: python
+
+    >>> hdulist = galex.serialize()
+    >>> hdulist.info()
+    Filename: (No file associated with this HDUList)
+    No.    Name      Ver    Type      Cards   Dimensions   Format
+    0  PRIMARY       1 PrimaryHDU       4   ()      
+    1                1 BinTableHDU     15   71002R x 1C   ['1J']   
+
+The hdulist contains two tables, a primary and a binary one. The MOC is stored in the second one.
+It consists of a list of the uniq numbers representing the set of HEALPix cells contained in the MOC (hence the 1d table).
+
+If you want to write it, just call the `~regions.MOCSkyRegion.write` method with the path
+you wish to save the MOC instance to.

docs/moc/plot_contains.py

@@ -0,0 +1,47 @@
+import numpy as np
+
+from astropy.utils.data import get_pkg_data_filename
+from regions import MOCSkyRegion, WCS
+from astropy.coordinates import Angle, SkyCoord
+import astropy.units as u
+# Load a MOC, e.g. the MOC of GALEXGR6-AIS-FUV
+filename = get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions')
+moc = MOCSkyRegion.from_fits(filename)
+
+# Generate random sky coordinates
+ra = np.random.randint(low=0, high=360, size=1000) * u.deg
+dec = np.random.randint(low=-90, high=90, size=1000) * u.deg
+# Get the mask of the sky coordinates contained in the MOCSkyRegion instance
+inside_mask = moc.contains(ra, dec)
+
+coords_inside = SkyCoord(ra[inside_mask], dec[inside_mask])
+coords_outside = SkyCoord(ra[~inside_mask], dec[~inside_mask])
+
+# Plot the MOC using matplotlib
+import matplotlib.pyplot as plt
+fig = plt.figure(111, figsize=(10, 10))
+# Define a WCS as a context
+with WCS(fig, 
+        fov=100 * 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)
+    # Cast the MOC sky region to a MOC pixel region
+    reg = moc.to_pixel(wcs)
+    # Call the plot method of MOC pixel region
+    reg.plot(ax=ax, alpha=0.5, fill=True, linewidth=0, color='r')
+    # Project the sky coordinates to the image system
+    from astropy.wcs.utils import skycoord_to_pixel
+    x_in, y_in = skycoord_to_pixel(coords_inside, wcs=wcs)
+    x_out, y_out = skycoord_to_pixel(coords_outside, wcs=wcs)
+    plt.scatter(x_out, y_out, s=16, c='black', alpha=0.5, marker='^', zorder=2, label='outside')
+    plt.scatter(x_in, y_in, s=16, c='green', alpha=0.5, marker='^', zorder=3, label='inside')
+
+plt.legend()
+plt.title("MOCSkyRegion of GALEX")
+plt.xlabel('ra')
+plt.ylabel('dec')
+plt.grid(color="black", linestyle="dotted")
+plt.show()

docs/moc/plot_logical_operations.py

@@ -0,0 +1,33 @@
+from astropy.utils.data import get_pkg_data_filename
+from regions import MOCSkyRegion, WCS
+from astropy.coordinates import Angle, SkyCoord
+import astropy.units as u
+# Load a MOC, e.g. the MOC of GALEXGR6-AIS-FUV
+galex = MOCSkyRegion.from_fits(get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions'))
+sdss9 = MOCSkyRegion.from_fits(get_pkg_data_filename('shapes/tests/data/P-SDSS9-r.fits', package='regions'))
+
+# Compute the union and intersection
+intersection = galex.intersection(sdss9)
+union = galex.union(sdss9)
+# Plot the MOC using matplotlib
+import matplotlib.pyplot as plt
+fig = plt.figure(111, figsize=(10, 10))
+# Define a WCS as a context
+with WCS(fig, 
+        fov=100 * 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)
+    union.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, linewidth=0, color='r', label="Union")
+    union.border(ax=ax, wcs=wcs, color='r')
+    intersection.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, linewidth=0, color='g', label="Intersection")
+    intersection.border(ax=ax, wcs=wcs, color='g')
+    ax.legend()
+
+plt.title("Logical operations between GALEX and SDSS9")
+plt.xlabel('ra')
+plt.ylabel('dec')
+plt.grid(color="black", linestyle="dotted")
+plt.show()

docs/moc/plot_moc.py

@@ -0,0 +1,32 @@
+# Example illustrating how to:
+# - Load a MOC from a FITS file
+# - Plot a MOC defining a WCS
+
+from astropy.utils.data import get_pkg_data_filename
+from regions import MOCSkyRegion, WCS
+from astropy.coordinates import Angle, SkyCoord
+import astropy.units as u
+# Load a MOC, e.g. the MOC of GALEXGR6-AIS-FUV
+filename = get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions')
+moc = MOCSkyRegion.from_fits(filename)
+# Plot the MOC using matplotlib
+import matplotlib.pyplot as plt
+fig = plt.figure(111, figsize=(10, 10))
+# Define a WCS as a context
+with WCS(fig, 
+        fov=100 * 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)
+    # Cast the MOC sky region to a MOC pixel region
+    reg = moc.to_pixel(wcs)
+    # Call the plot method of MOC pixel region
+    reg.plot(ax=ax, alpha=0.5, fill=True, linewidth=0.5, color='r')
+
+plt.title("MOCSkyRegion of GALEX")
+plt.xlabel('ra')
+plt.ylabel('dec')
+plt.grid(color="black", linestyle="dotted")
+plt.show()