diff --git a/.circleci/config.yml b/.circleci/config.yml
index e52bd359..dafb84c0 100644
--- a/.circleci/config.yml
+++ b/.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
diff --git a/.gitignore b/.gitignore
index 7aa981ba..fc94835d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,6 +28,9 @@ docs/_build
# Pycharm editor project files
.idea
+# VSCode editor project files
+.vscode/
+
# Packages/installer info
*.egg
*.egg-info
diff --git a/.rtd-environment.yml b/.rtd-environment.yml
index 97288507..249ca486 100644
--- a/.rtd-environment.yml
+++ b/.rtd-environment.yml
@@ -9,6 +9,7 @@ dependencies:
- numpy
- cython
- matplotlib
+ - spherical-geometry
- sphinx
- astropy
- shapely
diff --git a/.travis.yml b/.travis.yml
index 24cc92a2..60a4d90b 100644
--- a/.travis.yml
+++ b/.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
diff --git a/CHANGES.rst b/CHANGES.rst
index 237a8e9c..49ebb4d4 100644
--- a/CHANGES.rst
+++ b/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 `__. 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)
================
diff --git a/appveyor.yml b/appveyor.yml
index 3ad5eee2..50a21f3b 100644
--- a/appveyor.yml
+++ b/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
diff --git a/docs/index.rst b/docs/index.rst
index dafbe28c..7b57c6a8 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -32,6 +32,7 @@ User Documentation
installation
getting_started
shapes
+ moc/moc
contains
compound
masks
diff --git a/docs/moc/moc.rst b/docs/moc/moc.rst
new file mode 100644
index 00000000..48fb2a35
--- /dev/null
+++ b/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 `__ 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] `
+- 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 a cone (will be implemented in a next version).
+
+2. Perform fast logical operations between `~regions.MOCSkyRegion` objects: :ref:`[4] `
+
+- 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] `
+
+More informations about MOCs can be found by reading the following `IVOA paper `_.
+
+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] ` 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.
diff --git a/docs/moc/plot_contains.py b/docs/moc/plot_contains.py
new file mode 100644
index 00000000..c7397b01
--- /dev/null
+++ b/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()
\ No newline at end of file
diff --git a/docs/moc/plot_logical_operations.py b/docs/moc/plot_logical_operations.py
new file mode 100644
index 00000000..a3f88fbd
--- /dev/null
+++ b/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()
\ No newline at end of file
diff --git a/docs/moc/plot_moc.py b/docs/moc/plot_moc.py
new file mode 100644
index 00000000..15e21722
--- /dev/null
+++ b/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()
\ No newline at end of file
diff --git a/docs/moc/plot_moc_no_cast.py b/docs/moc/plot_moc_no_cast.py
new file mode 100644
index 00000000..6c2f3fa2
--- /dev/null
+++ b/docs/moc/plot_moc_no_cast.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)
+ # Add the MOC cells to the axe
+ moc.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, linewidth=0, color='r')
+ # Add the perimeter of the MOC to the axe
+ moc.border(ax=ax, wcs=wcs, alpha=0.5, fill=True, color='r')
+
+plt.title("MOCSkyRegion of GALEX")
+plt.xlabel('ra')
+plt.ylabel('dec')
+plt.grid(color="black", linestyle="dotted")
+plt.show()
\ No newline at end of file
diff --git a/docs/moc/plot_polygon.py b/docs/moc/plot_polygon.py
new file mode 100644
index 00000000..67764447
--- /dev/null
+++ b/docs/moc/plot_polygon.py
@@ -0,0 +1,52 @@
+from regions import MOCSkyRegion, WCS
+from astropy.coordinates import Angle, SkyCoord
+import astropy.units as u
+
+import numpy as np
+# The set of points delimiting the polygon in deg
+vertices = np.array([[18.36490956, 5. ],
+ [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. ]])
+skycoord = SkyCoord(vertices, unit="deg", frame="icrs")
+# A sky coordinate we know that belongs to the inside of the MOC
+inside = SkyCoord(ra=10, dec=5, unit="deg", frame="icrs")
+moc = MOCSkyRegion.from_polygon_skycoord(skycoord, max_depth=9, inside=inside)
+
+# Plot the MOC using matplotlib
+import matplotlib.pyplot as plt
+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('MOCSkyRegion from a polygon')
+plt.grid(color="black", linestyle="dotted")
+plt.show()
diff --git a/docs/shapes.rst b/docs/shapes.rst
index 6bc77d01..5fc26a7c 100644
--- a/docs/shapes.rst
+++ b/docs/shapes.rst
@@ -174,6 +174,22 @@ currently supported.
.. _sh-wcs:
+* `~regions.MOCSkyRegion`
+
+.. code-block:: python
+
+ >>> from astropy.coordinates import SkyCoord
+ >>> from regions import MOCSkyRegion
+ >>> # From an astropy skycoord
+ ... skycoord = SkyCoord([1, 2, 2], [1, 1, 2], unit='deg')
+ ... moc = MOCSkyRegion.from_skycoord(skycoord=skycoord, max_level=5)
+ >>> # From astropy lon and lat quantities
+ ... moc = MOCSkyRegion.from_lonlat(lon=skycoord.icrs.ra, lat=skycoord.icrs.dec, max_level=5)
+ >>> # From a FITS file
+ ... moc = MOCSkyRegion.from_fits()
+ >>> # From a python dict storing the HEALPix cell indexes
+ >>> moc = MOCSkyRegion.from_json(data={'2': [0, 6], '5': [1, 3, 5, ...], ...})
+
Transformations
---------------
diff --git a/regions/_geometry/pnpoly.pyx b/regions/_geometry/pnpoly.pyx
index b1c248ef..3fc58458 100644
--- a/regions/_geometry/pnpoly.pyx
+++ b/regions/_geometry/pnpoly.pyx
@@ -69,6 +69,36 @@ cimport numpy as np
DTYPE_BOOL = np.bool
ctypedef np.uint8_t DTYPE_BOOL_t
+def points_in_polygons(np.ndarray[DTYPE_t, ndim=1] x,
+ np.ndarray[DTYPE_t, ndim=1] y,
+ np.ndarray[DTYPE_t, ndim=2] vx,
+ np.ndarray[DTYPE_t, ndim=2] vy):
+ cdef int i, n
+ cdef np.ndarray[np.uint8_t, ndim=1, cast=True] result
+
+ n = x.shape[0]
+
+ result = np.zeros(n, DTYPE_BOOL)
+
+ for i in range(n):
+ result[i] = point_in_polygons(x[i], y[i], vx, vy)
+
+ return result
+
+cdef int point_in_polygons(double x, double y,
+ np.ndarray[DTYPE_t, ndim=2] vx,
+ np.ndarray[DTYPE_t, ndim=2] vy):
+ cdef int i, j, num_poly
+ cdef np.uint8_t in_poly
+
+ num_poly = vx.shape[0]
+
+ for i in range(num_poly):
+ in_poly = point_in_polygon(x, y, vx[i], vy[i])
+ if in_poly:
+ return 1
+
+ return 0
def points_in_polygon(np.ndarray[DTYPE_t, ndim=1] x,
np.ndarray[DTYPE_t, ndim=1] y,
diff --git a/regions/shapes/__init__.py b/regions/shapes/__init__.py
index 28f24a4e..14a79165 100644
--- a/regions/shapes/__init__.py
+++ b/regions/shapes/__init__.py
@@ -9,3 +9,4 @@
from .annulus import *
from .line import *
from .text import *
+from .moc import *
diff --git a/regions/shapes/moc/__init__.py b/regions/shapes/moc/__init__.py
new file mode 100644
index 00000000..da505eb8
--- /dev/null
+++ b/regions/shapes/moc/__init__.py
@@ -0,0 +1,8 @@
+from .core import MOCSkyRegion, MOCPixelRegion
+from .plot.wcs import WCS
+
+__all__ = [
+ 'MOCSkyRegion',
+ 'MOCPixelRegion',
+ 'WCS'
+]
diff --git a/regions/shapes/moc/boundaries.py b/regions/shapes/moc/boundaries.py
new file mode 100644
index 00000000..93c8854a
--- /dev/null
+++ b/regions/shapes/moc/boundaries.py
@@ -0,0 +1,189 @@
+import numpy as np
+# A python module handling graph manipulations
+import networkx as nx
+
+from astropy_healpix import HEALPix
+
+from astropy.coordinates import ICRS, SkyCoord
+from astropy.wcs.utils import skycoord_to_pixel
+
+from astropy_healpix import level_to_nside, nside_to_npix
+
+class Boundaries():
+ @staticmethod
+ def get(moc, depth):
+ boundaries_l = []
+
+ # Get the ipixels of the MOC at the deepest depth
+ hp, ipixels = Boundaries._compute_HEALPix_indices(moc, depth)
+
+ # Split the global MOC graph into all its non connected subgraphs.
+ graph_boundaries = Boundaries._compute_graph_HEALPix_boundaries(hp, ipixels)
+ boundaries_l.extend(Boundaries._retrieve_skycoords(graph_boundaries))
+
+ return boundaries_l
+
+ @staticmethod
+ def _compute_HEALPix_indices(m, depth):
+ moc = m
+ if depth:
+ if m.max_depth > depth:
+ moc = m.degrade_to_depth(depth)
+
+ max_depth = moc.max_depth
+ nside = level_to_nside(max_depth)
+ hp = HEALPix(nside=nside, order='nested', frame=ICRS())
+ ipixels = moc._best_res_pixels()
+
+ # Take the complement if the MOC covers more than half of the sky => the perimeter(MOC) = perimeter(complement(MOC))
+ # but we process less hpx cells
+ num_ipixels = nside_to_npix(nside)
+ sky_fraction = ipixels.shape[0] / float(num_ipixels)
+
+ #if sky_fraction > 0.5:
+ # ipixels_all = np.arange(num_ipixels)
+ # ipixels = np.setdiff1d(ipixels_all, ipixels, assume_unique=True)
+
+ return hp, ipixels
+
+ @staticmethod
+ def _compute_graph_HEALPix_boundaries(hp, ipixels):
+ def insert_edge(G, l1, l2, p1_lon, p1_lat, p2_lon, p2_lat):
+ # Nodes are indexed by str(skycoord). When getting ordered nodes, one can retrieve back the skycoord instance
+ # by accessing the python dict `pts_d`.
+ try:
+ # Avoid the special case where holes are touching to each other
+ # 'x' belongs to the MOC
+ # ' ' is part of the holes in the MOC
+ # |xxx
+ # |xxx
+ # ---A---
+ # xxx|
+ # xxx|
+ #
+ # If this case occurs we split the node A into 2. One is attached to the bottom left graph and the other to the
+ # top right one. When computing the MST (minimal spanning tree) from a graph, we need our graphs to have
+ # only nodes of degrees 1 or 2 (i.e. to be lines).
+ if G.degree[l1] >= 2:
+ l1 += '_'
+ except:
+ pass
+
+ try:
+ if G.degree[l2] >= 2:
+ l2 += '_'
+ except:
+ pass
+ # Set the skycoord instance as an attribute of the nodes
+ G.add_node(l1, ra=p1_lon, dec=p1_lat)
+ G.add_node(l2, ra=p2_lon, dec=p2_lat)
+ G.add_edge(l1, l2)
+
+ # Phase 1: Retrieve the ipixels located at the border of
+ # this connexe MOC component
+ ipix_neigh = hp.neighbours(ipixels)[[0, 2, 4, 6], :]
+
+ r1 = np.in1d(ipix_neigh[0, :], ipixels)
+ r2 = np.in1d(ipix_neigh[1, :], ipixels)
+ r3 = np.in1d(ipix_neigh[2, :], ipixels)
+ r4 = np.in1d(ipix_neigh[3, :], ipixels)
+
+ isin = np.vstack((r1, r2, r3, r4))
+
+ num_neigh = isin.sum(axis=0)
+ border = num_neigh < 4
+
+ ipixels_border = ipixels[border]
+ isin_border = isin[:, border]
+
+ # Phase 2: Build the graph from the positions of the ipixels boundaries
+ ipix_lon, ipix_lat = hp.boundaries_lonlat(ipixels_border, step=1)
+
+ ipix_lon_deg = ipix_lon.deg
+ ipix_lat_deg = ipix_lat.deg
+
+ ipix_lon_repr = \
+ np.around(np.asarray(ipix_lon.reshape((1, -1))[0]), decimals=6).tolist()
+ ipix_lat_repr = \
+ np.around(np.asarray(ipix_lat.reshape((1, -1))[0]), decimals=6).tolist()
+
+ west_border = ~isin_border[0, :]
+ south_border = ~isin_border[1, :]
+ east_border = ~isin_border[2, :]
+ north_border = ~isin_border[3, :]
+
+ E = nx.Graph()
+
+ for i in range(ipixels_border.shape[0]):
+ lon_deg = ipix_lon_deg[i]
+ lat_deg = ipix_lat_deg[i]
+
+ p0_lon = lon_deg[0]
+ p1_lon = lon_deg[1]
+ p2_lon = lon_deg[2]
+ p3_lon = lon_deg[3]
+
+ p0_lat = lat_deg[0]
+ p1_lat = lat_deg[1]
+ p2_lat = lat_deg[2]
+ p3_lat = lat_deg[3]
+
+ off = 4*i
+ off2 = 4*(i + 1)
+ repr_lon = ipix_lon_repr[off:off2]
+ repr_lat = ipix_lat_repr[off:off2]
+
+ s0 = str(repr_lon[0]) + ' ' + str(repr_lat[0])
+ s1 = str(repr_lon[1]) + ' ' + str(repr_lat[1])
+ s2 = str(repr_lon[2]) + ' ' + str(repr_lat[2])
+ s3 = str(repr_lon[3]) + ' ' + str(repr_lat[3])
+
+ # WEST border
+ if west_border[i]:
+ insert_edge(E, s1, s2, p1_lon, p1_lat, p2_lon, p2_lat)
+
+ # NORTH border
+ if north_border[i]:
+ insert_edge(E, s2, s3, p2_lon, p2_lat, p3_lon, p3_lat)
+
+ # EAST border
+ if east_border[i]:
+ insert_edge(E, s3, s0, p3_lon, p3_lat, p0_lon, p0_lat)
+
+ # SOUTH border
+ if south_border[i]:
+ insert_edge(E, s0, s1, p0_lon, p0_lat, p1_lon, p1_lat)
+
+ return E
+
+ @staticmethod
+ def _retrieve_skycoords(V):
+ coords_l = []
+ # Accessing the borders one by one. At this step, V_subgraphs contains a list of cycles
+ # (i.e. one describing the external border of the MOC component and several describing the holes
+ # found in the MOC component).
+ V_subgraphs = nx.connected_component_subgraphs(V)
+ for v in V_subgraphs:
+ # Compute the MST for each cycle
+ v = nx.convert_node_labels_to_integers(v)
+ mst = nx.minimum_spanning_tree(v)
+ # Get one end of the span tree by looping over its node and checking if the degree is one
+ src = None
+ for (node, deg) in mst.degree():
+ if deg == 1:
+ src = node
+ break
+
+ # Get the unordered lon and lat
+ ra = np.asarray(list(nx.get_node_attributes(v, 'ra').values()))
+ dec = np.asarray(list(nx.get_node_attributes(v, 'dec').values()))
+ coords = np.vstack((ra, dec)).T
+ # Get the ordering from the MST
+ ordering = np.asarray(list(nx.dfs_preorder_nodes(mst, src)))
+ # Order the coords
+ coords = coords[ordering]
+ # Get a skycoord containing N coordinates computed from the Nx2 `coords` array
+ coords = SkyCoord(coords, unit="deg")
+ coords_l.append(coords)
+
+ return coords_l
diff --git a/regions/shapes/moc/core.py b/regions/shapes/moc/core.py
new file mode 100644
index 00000000..d19c9862
--- /dev/null
+++ b/regions/shapes/moc/core.py
@@ -0,0 +1,1231 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function
+from .py23_compat import range, int
+
+import numpy as np
+from astropy import units as u
+
+from astropy import wcs
+from astropy.io import fits
+from astropy.coordinates import ICRS, SkyCoord
+from astropy.table import Table
+from astropy.wcs.utils import pixel_to_skycoord, \
+ skycoord_to_pixel
+
+from astropy_healpix import HEALPix, \
+ nside_to_npix, \
+ level_to_nside, \
+ uniq_to_level_ipix
+
+from .interval_set import IntervalSet
+from .utils import trailing_zeros
+
+from .boundaries import Boundaries
+from .plot import fill
+from .plot import border
+from .plot import axis_viewport
+from .plot import culling_backfacing_cells
+
+from ...core import PixCoord, \
+ PixelRegion, \
+ SkyRegion, \
+ BoundingBox, \
+ RegionMask
+from ...core.attributes import RegionMeta, \
+ RegionVisual
+
+from ..._geometry.pnpoly import points_in_polygons
+from ..._geometry import polygonal_overlap_grid
+
+__all__ = [
+ 'MOCSkyRegion',
+ 'MOCPixelRegion',
+]
+
+class MOCPixelRegion(PixelRegion):
+ """
+ A MOC (Multi-Order Coverage) map in pixel coordinates.
+
+ Parameters
+ ----------
+ sky_region : `~regions.MOCSkyRegion`
+ The `~regions.MOCSkyRegion` instance.
+ meta : `~regions.RegionMeta`, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual`, optional
+ A dictionary which stores the visual meta attributes of this region.
+ """
+
+ def __init__(self, wcs, sky_region, meta=None, visual=None):
+ self.moc = sky_region
+ self.meta = meta or RegionMeta()
+ self.visual = visual or RegionVisual()
+ self._repr_params = ('_interval_set',)
+ self.wcs = wcs
+ self.vertices = self._compute_vertices()
+
+ def _compute_vertices(self):
+ # Simplify the MOC for plotting purposes:
+ # 1. Degrade the MOC if the FOV is enough big so that we cannot see the smallest HEALPix cells.
+ # 2. For small FOVs, plot the MOC & POLYGONAL_MOC_FROM_FOV
+ moc = fill.build_plotting_moc(moc=self.moc, wcs=self.wcs)
+ # If the FOV contains no cells, then moc_to_plot (i.e. the intersection between the moc
+ # and the MOC created from the FOV polygon) will be empty.
+ # If it is the case, we exit the method without doing anything.
+ vertices = np.array([], dtype=float)
+ if moc.empty():
+ return vertices
+
+ depth_ipix_d = moc.serialize(format="json")
+ depth_ipix_clean_d = culling_backfacing_cells.from_moc(depth_ipix_d=depth_ipix_d, wcs=self.wcs)
+
+ for depth, ipix in depth_ipix_clean_d.items():
+ step = 1
+ depth = int(depth)
+ if depth < 3:
+ step = 2
+
+ nside = level_to_nside(depth)
+ hp = HEALPix(order="nested", nside=nside, frame=ICRS())
+ ipix_boundaries = hp.boundaries_skycoord(ipix, step=step)
+ # Projection on the given WCS
+ xp, yp = skycoord_to_pixel(ipix_boundaries, wcs=self.wcs)
+
+ patches = np.vstack((xp.flatten(), yp.flatten())).T
+ patches = patches.reshape((-1, step*4, 2))
+
+ if vertices.size == 0:
+ vertices = patches
+ else:
+ vertices = np.append(vertices, patches, axis=0)
+
+ return vertices
+
+ def as_artist(self, origin=(0, 0), **kw_mpl_pathpatch):
+ """
+ Matplotlib patch object for this region (`matplotlib.patches.PathPatch`)
+
+ Parameters:
+ -----------
+ origin : `~numpy.ndarray`, optional
+ The ``(x, y)`` pixel position of the origin of the displayed image.
+ Default is (0, 0).
+ kwargs: dict
+ All keywords that a `~matplotlib.patches.PathPatch` object accepts
+
+ Returns
+ -------
+ patch : `~matplotlib.patches.PathPatch`
+ Matplotlib path patch
+
+ Examples
+ --------
+ >>> 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")
+ """
+ pathpatch_mpl = fill.fill(moc=self.moc, wcs=self.wcs, origin=origin, **kw_mpl_pathpatch)
+ return pathpatch_mpl
+
+ def plot(self, ax=None, origin=(0, 0), **kwargs):
+ """
+ Calls ``as_artist`` method forwarding all kwargs and adds patch
+ to given axis.
+
+ Parameters
+ ----------
+ origin : array_like, optional
+ The ``(x, y)`` pixel position of the origin of the displayed image.
+ Default is (0, 0).
+ ax : `~matplotlib.axes.Axes`, optional
+ Axis
+ perimeter : bool, optional
+ Plot the perimeter of the MOCPixelRegion. Default to True.
+ kwargs : `dict`
+ keywords that a `~matplotlib.patches.Patch` accepts
+
+ Returns
+ -------
+ ax : `~matplotlib.axes.Axes`
+ Axes on which the patch is added.
+ """
+ PixelRegion.plot(self, ax=ax, origin=origin, **kwargs)
+ axis_viewport.set(ax=ax, wcs=self.wcs)
+ return ax
+
+ def to_sky(self, wcs):
+ return self.moc
+
+ @property
+ def bounding_box(self):
+ xmin = np.min(self.vertices[:, :, 0])
+ xmax = np.max(self.vertices[:, :, 0])
+ ymin = np.min(self.vertices[:, :, 1])
+ ymax = np.max(self.vertices[:, :, 1])
+ return BoundingBox.from_float(xmin, xmax, ymin, ymax)
+
+ @property
+ def area(self):
+ raise NotImplementedError
+
+ def to_mask(self, mode='center', subpixels=5):
+ self._validate_mode(mode, subpixels)
+
+ if mode == 'center':
+ mode = 'subpixels'
+ subpixels = 1
+
+ if mode == 'subpixels':
+ use_exact = 0
+ else:
+ use_exact = 1
+
+ # Find bounding box and mask size
+ bbox = self.bounding_box
+ ny, nx = bbox.shape
+
+ # Find position of pixel edges and recenter so that circle is at origin
+ xmin = float(bbox.ixmin) - 0.5
+ xmax = float(bbox.ixmax) - 0.5
+ ymin = float(bbox.iymin) - 0.5
+ ymax = float(bbox.iymax) - 0.5
+
+ vx, vy = self.vertices[:, :, 0], self.vertices[:, :, 1]
+
+ # Loop over all the projeted HEALPix cells to get their overlap grids
+ buffer = np.zeros(shape=(ny, nx))
+ for i in range(vx.shape[0]):
+ ovlp_grid = polygonal_overlap_grid(
+ xmin, xmax, ymin, ymax,
+ nx, ny, vx[i], vy[i],
+ use_exact, subpixels,
+ )
+ # Sum this overlap grid to the buffer
+ buffer += ovlp_grid
+
+ # Clip its values to the interval (0, 1)
+ mask = np.clip(a=buffer, a_min=0, a_max=1)
+ return RegionMask(mask, bbox=bbox)
+
+ def contains(self, pixcoord):
+ """
+ Checks whether a position or positions fall inside the region.
+
+ Parameters
+ ----------
+ pixcoord : `~regions.PixCoord`
+ The position or positions to check.
+ """
+ pixcoord = PixCoord._validate(pixcoord, 'pixcoord')
+ x = np.atleast_1d(np.asarray(pixcoord.x, dtype=float))
+ y = np.atleast_1d(np.asarray(pixcoord.y, dtype=float))
+
+ shape = x.shape
+ vx, vy = self.vertices[:, :, 0], self.vertices[:, :, 1]
+ mask = points_in_polygons(x.flatten(), y.flatten(), vx, vy).astype(bool)
+ in_poly = mask.reshape(shape)
+ if self.meta.get('include', True):
+ return in_poly
+ else:
+ return np.logical_not(in_poly)
+
+class MOCSkyRegion(SkyRegion):
+ """
+ Multi-order spatial coverage class
+
+ 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 `__ 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 `from_fits`).
+ - Directly from a list of HEALPix cells expressed either as a numpy structural array (see `from_cells`) or a simple
+ python dictionnary (see `from_json`).
+ - From a list of sky coordinates (see `from_skycoord`, `from_lonlat`).
+ - From a convex/concave polygon (see `from_polygon`).
+ - From a cone (will be implemented in a next version).
+
+ 2. Perform fast logical operations between `regions.MOCSkyRegion` objects:
+
+ - The `intersection`
+ - The `union`
+ - The `difference`
+ - The `complement`
+
+ 3. Plot the `regions.MOCSkyRegion` objects:
+
+ - Draw the MOC with its HEALPix cells (see `fill`)
+ - Draw the perimeter of a MOC (see `border`)
+
+ 4. Get the sky coordinates defining the border(s) of `regions.MOCSkyRegion` objects (see `get_boundaries`).
+
+ 5. Serialize `regions.MOCSkyRegion` objects to `astropy.io.fits.HDUList` or JSON dictionary and save it to a file.
+
+ Parameters
+ ----------
+ itv_s : `~regions.IntervalSet` object, optional
+ A N rows by 2 columns `~numpy.ndarray` storing the ranges of HEALPix cells contained in the MOC.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+ """
+ HPY_MAX_DEPTH = 29
+
+ def __init__(self, itv=None, meta=None, visual=None):
+ interval = IntervalSet() if itv is None else itv
+ self._itv = interval
+ self.meta = meta or RegionMeta()
+ self.visual = visual or RegionVisual()
+
+ def to_pixel(self, wcs):
+ """
+ Convert the `~regions.MOCSkyRegion` instance to a `~regions.MOCPixelRegion`
+
+ The MOC is projected into the image coordinates system using an `~astropy.wcs.WCS` instance.
+ The HEALPix cells that backface the projection are removed.
+
+ Parameters
+ ----------
+ wcs : `~astropy.wcs.WCS`
+ The world coordinate system <-> Image coordinate system projection
+
+ Returns
+ -------
+ pixel_region : `~regions.MOCPixelRegion`
+ The projection of the MOCSkyRegion instance.
+ """
+ return MOCPixelRegion(wcs=wcs, sky_region=self, meta=self.meta, visual=self.visual)
+
+ def __repr__(self):
+ return self._itv.__repr__()
+
+ def __eq__(self, other):
+ """
+ Test equality between two `~regions.MOCSkyRegion`
+
+ Parameters
+ ----------
+ other : `~regions.MOCSkyRegion`
+ The `~regions.MOCSkyRegion` instance to test the equality with.
+
+ Returns
+ -------
+ result : bool
+ True is the two `~regions.MOCSkyRegion` are equal.
+ """
+ if not isinstance(other, MOCSkyRegion):
+ raise TypeError('Cannot compare a MOCSkyRegion with a {0}'.format(type(other)))
+
+ return self._itv == other._itv
+
+ def empty(self):
+ """
+ Checks whether the `~regions.MOCSkyRegion` instance is empty
+
+ A MOC is empty when it contains no ranges of HEALPix cell indexes.
+
+ Returns
+ -------
+ result: bool
+ True if the MOC instance is empty.
+ """
+ return self._itv.empty()
+
+ @property
+ def max_depth(self):
+ """
+ Depth of the smallest HEALPix cells found in the MOC instance.
+ """
+ combo = int(0)
+ for iv in self._itv._data:
+ combo |= iv[0] | iv[1]
+
+ ret = MOCSkyRegion.HPY_MAX_DEPTH - trailing_zeros(combo)//2
+ if ret < 0:
+ ret = 0
+
+ return ret
+
+ def intersection(self, other, *args):
+ """
+ Intersection between the MOC instance and other MOCs.
+
+ Parameters
+ ----------
+ other : `regions.MOCSkyRegion`
+ The MOC used for performing the intersection with self.
+ args : `regions.MOCSkyRegion`
+ Other additional MOCs to perform the intersection with.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC.
+ """
+ itv = self._itv.intersection(other._itv)
+ for moc in args:
+ itv = itv.intersection(moc._itv)
+
+ return self.__class__(itv)
+
+ def union(self, other, *args):
+ """
+ Union between the MOC instance and other MOCs.
+
+ Parameters
+ ----------
+ other : `regions.MOCSkyRegion`
+ The MOC used for performing the union with self.
+ args : `regions.MOCSkyRegion`
+ Other additional MOCs to perform the union with.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC.
+ """
+ itv = self._itv.union(other._itv)
+ for moc in args:
+ itv = itv.union(moc._itv)
+
+ return self.__class__(itv)
+
+ def difference(self, other, *args):
+ """
+ Difference between the MOC instance and other MOCs.
+
+ Parameters
+ ----------
+ other : `regions.MOCSkyRegion`
+ The MOC used that will be substracted to self.
+ args : `regions.MOCSkyRegion`
+ Other additional MOCs to perform the difference with.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC.
+ """
+ itv = self._itv.difference(other._itv)
+ for moc in args:
+ itv = itv.difference(moc._itv)
+
+ return self.__class__(itv)
+
+ def complement(self):
+ """
+ Returns the complement of the MOC instance.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC.
+ """
+ res = []
+ itvs_l = sorted(self._itv._data.tolist())
+
+ if itvs_l[0][0] > 0:
+ res.append((0, itvs_l[0][0]))
+
+ last = itvs_l[0][1]
+
+ for itv in itvs_l[1:]:
+ res.append((last, itv[0]))
+ last = itv[1]
+
+ max_npix = 3 << 60
+
+ if last < max_npix:
+ res.append((last, max_npix))
+
+ itvs = np.asarray(res, dtype=np.int64)
+ itv = IntervalSet(itvs)
+ return self.__class__(itv)
+
+
+ @classmethod
+ def from_cells(cls, cells, meta=None, visual=None):
+ """
+ Creates a MOC from a numpy array representing the HEALPix cells.
+
+ Parameters
+ ----------
+ cells : `numpy.ndarray`
+ Must be a numpy structured array (See https://docs.scipy.org/doc/numpy-1.15.0/user/basics.rec.html).
+ The structure of a cell contains 3 attributes:
+
+ - A `ipix` value being a np.uint64
+ - A `depth` value being a np.uint32
+ - A `fully_covered` flag bit stored in a np.uint8
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ moc : `regions.MOCSkyRegion`
+ The MOC.
+ """
+ shift = (MOCSkyRegion.HPY_MAX_DEPTH - cells["depth"]) << 1
+
+ p1 = cells["ipix"]
+ p2 = cells["ipix"] + 1
+
+ itvs = np.vstack((p1 << shift, p2 << shift)).T
+
+ return cls(IntervalSet(itvs), meta, visual)
+
+ @classmethod
+ def from_json(cls, data, meta=None, visual=None):
+ """
+ Creates a `~regions.MOCSkyRegion` instance from a dictionary of HEALPix cell arrays indexed by their depth.
+
+ Parameters
+ ----------
+ data : dict(str : [int])
+ A dictionary of HEALPix cell arrays indexed by their depth.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ moc : `regions.MOCSkyRegion`
+ The MOC instance.
+
+ Examples
+ --------
+ >>> from regions import MOCSkyRegion
+ >>> moc = MOCSkyRegion.from_json({
+ ... "0": [0, 3, 4],
+ ... "1": [5, 6],
+ ... })
+ """
+ itvs = np.array([], dtype=np.int64)
+ for order, pix_l in data.items():
+ if len(pix_l) == 0:
+ continue
+ pix = np.array(pix_l, dtype=np.int64)
+ p1 = pix
+ p2 = pix + 1
+ shift = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - int(order))
+
+ itv = np.vstack((p1 << shift, p2 << shift)).T
+ if itvs.size == 0:
+ itvs = itv
+ else:
+ itvs = np.vstack((itvs, itv))
+
+ itv_s = IntervalSet(itvs)
+ return cls(itv_s, meta, visual)
+
+ def _uniq_pixels_iterator(self):
+ """
+ Generator giving the NUNIQ HEALPix pixels of the MOC.
+
+ Returns
+ -------
+ uniq :
+ the NUNIQ HEALPix pixels iterator
+ """
+ itvs_uniq_l = IntervalSet.to_uniq_itv_s(self._itv)._data
+ for uniq_iv in itvs_uniq_l:
+ for uniq in np.arange(uniq_iv[0], uniq_iv[1], dtype=np.int64):
+ yield uniq
+
+ @classmethod
+ def from_fits(cls, filename, meta=None, visual=None):
+ """
+ Creates a `~regions.MOCSkyRegion` from a FITS file.
+
+ The specified FITS file must store the MOC (i.e. the list of HEALPix cells it contains) in a binary HDU table.
+
+ Parameters
+ ----------
+ filename : str
+ The path to the FITS file.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ moc : `~regions.MOCSkyRegion`
+ The MOC instance.
+ """
+ table = Table.read(filename)
+
+ itvs = np.vstack((table['UNIQ'], table['UNIQ']+1)).T
+
+ nuniq_itv = IntervalSet(itvs)
+ itv = IntervalSet.from_uniq_itv_s(nuniq_itv)
+ return cls(itv, meta, visual)
+
+ @staticmethod
+ def _to_json(uniq):
+ """
+ Serializes a MOC to the JSON format.
+
+ Parameters
+ ----------
+ uniq : `~numpy.ndarray`
+ The array of HEALPix cells representing the MOC to serialize.
+
+ Returns
+ -------
+ result : dict(str : [int])
+ A dictionary of HEALPix cell lists indexed by their depth.
+ """
+ result = {}
+
+ depth, ipix = uniq_to_level_ipix(uniq)
+ min_depth = np.min(depth[0])
+ max_depth = np.max(depth[-1])
+
+ for d in range(min_depth, max_depth+1):
+ pix_index = np.where(depth == d)[0]
+ if pix_index.size:
+ # There are pixels belonging to the current order
+ ipix_depth = ipix[pix_index]
+ result[str(d)] = ipix_depth.tolist()
+
+ return result
+
+ def _to_fits(self, uniq_pix, optional_kw_dict=None):
+ """
+ Serializes a MOC to the FITS format.
+
+ Parameters
+ ----------
+ uniq_pix : `numpy.ndarray`
+ The array of HEALPix cells representing the MOC to serialize.
+ optional_kw_dict : dict
+ Optional keywords arguments added to the FITS header.
+
+ Returns
+ -------
+ thdulist : `astropy.io.fits.HDUList`
+ The list of HDU tables.
+ """
+ depth = self.max_depth
+ if depth <= 13:
+ fits_format = '1J'
+ else:
+ fits_format = '1K'
+
+ tbhdu = fits.BinTableHDU.from_columns(
+ fits.ColDefs([
+ fits.Column(name='UNIQ', format=fits_format, array=uniq_pix)
+ ]))
+ tbhdu.header['PIXTYPE'] = 'HEALPIX'
+ tbhdu.header['ORDERING'] = 'NUNIQ'
+ tbhdu.header['COORDSYS'] = ('C', 'reference frame (C=ICRS)')
+ tbhdu.header['MOCORDER'] = depth
+ tbhdu.header['MOCTOOL'] = 'MOCPy'
+ if optional_kw_dict:
+ for key in optional_kw_dict:
+ tbhdu.header[key] = optional_kw_dict[key]
+
+ thdulist = fits.HDUList([fits.PrimaryHDU(), tbhdu])
+ return thdulist
+
+ def serialize(self, format='fits', optional_kw_dict=None):
+ """
+ Serializes the MOC into a specific format.
+
+ Possible formats are FITS and JSON.
+
+ Parameters
+ ----------
+ format : str
+ 'fits' by default. The other possible choice is 'json'.
+ optional_kw_dict : dict
+ Optional keywords arguments added to the FITS header. Only used if ``format`` equals to 'fits'.
+
+ Returns
+ -------
+ result : `astropy.io.fits.HDUList` or JSON dictionary
+ The result of the serialization.
+ """
+ formats = ('fits', 'json')
+ if format not in formats:
+ raise ValueError('format should be one of %s' % (str(formats)))
+
+ uniq_l = []
+ for uniq in self._uniq_pixels_iterator():
+ uniq_l.append(uniq)
+
+ uniq = np.array(uniq_l, dtype=np.int64)
+
+ if format == 'fits':
+ result = self._to_fits(uniq_pix=uniq,
+ optional_kw_dict=optional_kw_dict)
+ else:
+ # json format serialization
+ result = self.__class__._to_json(uniq=uniq)
+
+ return result
+
+ def write(self, path, format='fits', optional_kw_dict=None):
+ """
+ Writes the MOC to a file.
+
+ Format can be 'fits' or 'json', though only the fits format is officially supported by the IVOA.
+
+ Parameters
+ ----------
+ path : str, optional
+ The path to the file to save the MOC in.
+ format : str, optional
+ The format in which the MOC will be serialized before being saved. Possible formats are "fits" or "json".
+ By default, ``format`` is set to "fits".
+ optional_kw_dict : optional
+ Optional keywords arguments added to the FITS header. Only used if ``format`` equals to 'fits'.
+ """
+ serialization = self.serialize(format=format, optional_kw_dict=optional_kw_dict)
+ if format == 'fits':
+ serialization.writeto(path, overwrite=True)
+ else:
+ import json
+ with open(path, 'w') as h:
+ h.write(json.dumps(serialization, sort_keys=True, indent=2))
+
+ def degrade_to_depth(self, new_depth):
+ """
+ Degrades the MOC instance to a new, less precise, MOC.
+
+ The maximum depth (i.e. the depth of the smallest HEALPix cells that can be found in the MOC) of the
+ degraded MOC is set to ``new_depth``.
+
+ Parameters
+ ----------
+ new_depth : int
+
+ Returns
+ -------
+ moc : `regions.MOCSkyRegion`
+ The degraded MOC.
+ """
+ shift = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - new_depth)
+ ofs = (int(1) << shift) - 1
+ mask = ~ofs
+ adda = int(0)
+ addb = ofs
+ iv_set = []
+
+ for iv in self._itv._data:
+ a = (iv[0] + adda) & mask
+ b = (iv[1] + addb) & mask
+ if b > a:
+ iv_set.append((a, b))
+
+ return self.__class__(IntervalSet(np.asarray(iv_set, dtype=np.int64)))
+
+ def _best_res_pixels(self):
+ """
+ Returns a numpy array of all the HEALPix indexes contained in the MOC at its max depth.
+
+ Returns
+ -------
+ result : `~numpy.ndarray`
+ The array of HEALPix at ``max_depth``
+ """
+ factor = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - self.max_depth)
+ pix_l = []
+ for iv in self._itv._data:
+ for val in np.arange(iv[0] >> factor, iv[1] >> factor, dtype=np.int64):
+ pix_l.append(val)
+
+ return np.asarray(pix_l, dtype=np.int64)
+
+ def contains(self, ra, dec, keep_inside=True):
+ """
+ Returns a boolean mask array of the positions lying inside (or outside) the MOC instance.
+
+ Parameters
+ ----------
+ ra : `astropy.units.Quantity`
+ Right ascension array
+ dec : `astropy.units.Quantity`
+ Declination array
+
+ Returns
+ -------
+ array : `~numpy.ndarray`
+ A boolean numpy array telling which positions are inside the MOC depending on the
+ value of ``self.meta['include']``.
+ """
+ depth = self.max_depth
+
+ nside = level_to_nside(depth)
+ npix = nside_to_npix(nside)
+ mask = np.zeros(npix, dtype=bool)
+
+ ipix = self._best_res_pixels()
+ mask[ipix] = True
+
+ if self.meta.get('include', False):
+ mask = np.logical_not(mask)
+
+ # Retrieve the cells containing the sky positions
+ hp = HEALPix(nside=nside, order='nested')
+ ipix = hp.lonlat_to_healpix(ra, dec)
+
+ return mask[ipix]
+
+ def add_neighbours(self):
+ """
+ Extends the MOC instance so that it includes the HEALPix cells touching its border.
+
+ The depth of the HEALPix cells added at the border is equal to the maximum depth of the MOC instance.
+
+ Returns
+ -------
+ moc : `regions.MOCSkyRegion`
+ A new MOC instance that extends of one degree the initial MOC.
+ """
+ # Retrieve the ipixels at the deepest level
+ ipix = self._best_res_pixels()
+ max_depth = self.max_depth
+ nside = level_to_nside(max_depth)
+
+ # Compute the neighbours of the ipixels retrieved
+ hp = HEALPix(nside=nside, order='nested')
+ ipix_neigh = hp.neighbours(ipix)[[0, 2, 4, 6], :]
+ ipix_neigh = ipix_neigh.reshape((1, -1))[0]
+
+ # Get the union of the ipixels with their neighbours
+ res = np.union1d(ipix, ipix_neigh)
+
+ shift = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - max_depth)
+ itvs = np.vstack((res << shift, (res + 1) << shift)).T
+
+ itv_s = IntervalSet(itvs)
+ return MOCSkyRegion(itv_s)
+
+ def remove_neighbours(self):
+ """
+ Removes from the MOC instance the HEALPix cells located at its border.
+
+ The depth of the HEALPix cells removed is equal to the maximum depth of the MOC instance.
+
+ Returns
+ -------
+ moc : `regions.MOCSkyRegion`
+ A new MOC instance from which the HEALPix cells located at the border
+ of the initial MOC have been removed.
+ """
+ # Retrieve the ipixels at the deepest level
+ ipix = self._best_res_pixels()
+ max_depth = self.max_depth
+ nside = level_to_nside(max_depth)
+
+ # Retrieve the ipixels being at the border
+ hp = HEALPix(nside=nside, order='nested')
+ ipix_neigh = hp.neighbours(ipix)[[0, 2, 4, 6], :]
+
+ r1 = np.in1d(ipix_neigh[0, :], ipix)
+ r2 = np.in1d(ipix_neigh[1, :], ipix)
+ r3 = np.in1d(ipix_neigh[2, :], ipix)
+ r4 = np.in1d(ipix_neigh[3, :], ipix)
+
+ mask = np.vstack((r1, r2, r3, r4))
+
+ num_neigh = mask.sum(axis=0)
+ border = num_neigh < 4
+
+ # Get the ipixels which are not at the border
+ res = ipix[~border]
+
+ shift = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - max_depth)
+ itvs = np.vstack((res << shift, (res + 1) << shift)).T
+
+ itv_s = IntervalSet(itvs)
+ return MOCSkyRegion(itv_s)
+
+ def fill(self, ax, wcs, **kw_mpl_pathpatch):
+ """
+ Draws the MOC on a matplotlib axis.
+
+ This performs the projection of the cells from the world coordinate system to the pixel image coordinate system.
+ You are able to specify various styling kwargs for `matplotlib.patches.PathPatch`
+ (see the `list of valid keywords `__).
+
+ Parameters
+ ----------
+ ax : `matplotlib.axes.Axes`
+ Matplotlib axis.
+ wcs : `astropy.wcs.WCS`
+ WCS defining the World system <-> Image system projection.
+ kw_mpl_pathpatch
+ Plotting arguments for `matplotlib.patches.PathPatch`.
+
+ Examples
+ --------
+ >>> 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=(15, 15))
+ >>> # Define a WCS as a context
+ >>> with WCS(fig,
+ ... fov=50 * u.deg,
+ ... center=SkyCoord(0, 20, 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 giving the matplotlib axe and the `~astropy.wcs.WCS` object.
+ ... # We will set the matplotlib keyword linewidth to 0 so that it does not plot
+ ... # the border of each HEALPix cell.
+ ... # The color can also be specified along with an alpha value.
+ ... moc.fill(ax=ax, wcs=wcs, linewidth=0, alpha=0.5, fill=True, color="green")
+ >>> plt.xlabel('ra')
+ >>> plt.ylabel('dec')
+ >>> plt.grid(color="black", linestyle="dotted")
+ """
+ mpl_pathpatch = fill.fill(moc=self, wcs=wcs, origin=(0, 0), **kw_mpl_pathpatch)
+ ax.add_patch(mpl_pathpatch)
+ axis_viewport.set(ax=ax, wcs=wcs)
+
+ def border(self, ax, wcs, **kw_mpl_pathpatch):
+ """
+ Draws the MOC border(s) on a matplotlib axis.
+
+ This performs the projection of the sky coordinates defining the perimeter of the MOC to the pixel image coordinate system.
+ You are able to specify various styling kwargs for `matplotlib.patches.PathPatch`
+ (see the `list of valid keywords `__).
+
+ Parameters
+ ----------
+ ax : `matplotlib.axes.Axes`
+ Matplotlib axis.
+ wcs : `astropy.wcs.WCS`
+ WCS defining the World system <-> Image system projection.
+ kw_mpl_pathpatch
+ Plotting arguments for `matplotlib.patches.PathPatch`
+
+ Examples
+ --------
+ >>> 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=(15, 15))
+ >>> # Define a WCS as a context
+ >>> with WCS(fig,
+ ... fov=50 * u.deg,
+ ... center=SkyCoord(0, 20, 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 border giving the matplotlib axe and the `~astropy.wcs.WCS` object.
+ ... moc.border(ax=ax, wcs=wcs, alpha=0.5, color="red")
+ >>> plt.xlabel('ra')
+ >>> plt.ylabel('dec')
+ >>> plt.grid(color="black", linestyle="dotted")
+ """
+ border.border(self, ax, wcs, **kw_mpl_pathpatch)
+
+ def get_boundaries(self, depth=None):
+ """
+ Returns the sky coordinates defining the border(s) of the MOC.
+
+ The border(s) are expressed as a list of SkyCoord.
+ Each SkyCoord refers to the coordinates of one border of the MOC (i.e.
+ either a border of a connexe MOC part or a border of a hole
+ located in a connexe MOC part).
+
+ Parameters
+ ----------
+ depth : int
+ The depth of the MOC before computing its boundaries.
+ A shallow depth leads to a faster computation.
+ By default the maximum depth of the MOC is taken.
+
+ Return
+ ------
+ coords: [`~astropy.coordinates.SkyCoord`]
+ A list of `~astropy.coordinates.SkyCoord` each describing one border.
+ """
+ return Boundaries.get(self, depth)
+
+ @classmethod
+ def from_image(cls, header, max_depth, mask=None, meta=None, visual=None):
+ """
+ Creates a `~regions.MOCSkyRegion` instance from an image stored as a FITS file.
+
+ Parameters
+ ----------
+ header : `astropy.io.fits.Header`
+ FITS header containing all the infos of where the image is located (position, size, etc...)
+ max_depth : int
+ The moc resolution.
+ mask : `numpy.ndarray`, optional
+ A boolean array of the same size of the image where pixels having the value 1 are part of
+ the final MOC and pixels having the value 0 are not.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ moc : `~regions.MOCSkyRegion`
+ The MOC instance.
+ """
+ # Get the image dimensions
+ height = header['NAXIS2']
+ width = header['NAXIS1']
+
+ # Get the image WCS to retrieve back its world coordinates
+ w = wcs.WCS(header)
+
+ # Get the pixels coordinates
+ if mask is not None:
+ # Get the pixels in the mask if available
+ y, x = np.where(mask)
+ pix = np.dstack((x, y))[0]
+ else:
+ # Get a uniform grid of pixels whose boundaries are the dimensions of the image.
+ step = 1
+ x, y = np.mgrid[0.5:(width + 0.5 + step):step, 0.5:(height + 0.5 + step):step]
+ pix = np.dstack((x.ravel(), y.ravel()))[0]
+
+ # Retrieve back the world coordinates using the image WCS
+ world_pix = w.wcs_pix2world(pix, 1)
+
+ # Call the from_lonlat method to create the MOC from the world coordinates.
+ return MOCSkyRegion.from_lonlat(
+ lon=world_pix[:, 0] * u.deg,
+ lat=world_pix[:, 1] * u.deg,
+ max_depth=max_depth,
+ meta=meta,
+ visual=visual)
+
+ @classmethod
+ def from_fits_images(cls, path_l, max_depth, meta=None, visual=None):
+ """
+ Loads a MOC from a set of FITS file images.
+
+ Parameters
+ ----------
+ path_l : [str]
+ A list of path where the fits image are located.
+ max_depth : int
+ The MOC resolution.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ moc : `regions.MOCSkyRegion`
+ The union of all the MOCs created from the paths found in ``path_l``.
+ """
+ moc = MOCSkyRegion()
+ for path in path_l:
+ header = fits.getheader(path)
+ current_moc = MOCSkyRegion.from_image(header=header, max_depth=max_depth, meta=meta, visual=meta)
+ moc = moc.union(current_moc)
+
+ return moc
+
+ @classmethod
+ def from_skycoord(cls, skycoord, max_depth, meta=None, visual=None):
+ """
+ Creates a `~regions.MOCSkyRegion` from a `~astropy.coordinates.SkyCoord`.
+
+ Parameters
+ ----------
+ skycoord : `astropy.coordinates.SkyCoord`
+ The sky coordinates that will belong to the MOC.
+ max_depth : int
+ The depth of the smallest HEALPix cells contained in the MOC. Must be <= 29.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ moc : `~regions.MOCSkyRegion`
+ The MOC instance.
+ """
+ nside = level_to_nside(max_depth)
+ hp = HEALPix(nside=(1 << max_depth), order='nested')
+ ipix = hp.lonlat_to_healpix(skycoord.icrs.ra, skycoord.icrs.dec)
+
+ shift = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - max_depth)
+ itvs = np.vstack((ipix << shift, (ipix + 1) << shift)).T
+
+ itv = IntervalSet(itvs)
+ return cls(itv, meta, visual)
+
+ @classmethod
+ def from_lonlat(cls, lon, lat, max_depth, meta=None, visual=None):
+ """
+ Creates a MOC from astropy lon, lat `astropy.units.Quantity`.
+
+ Parameters
+ ----------
+ lon : `astropy.units.Quantity`
+ The longitudes of the sky coordinates belonging to the MOC.
+ lat : `astropy.units.Quantity`
+ The latitudes of the sky coordinates belonging to the MOC.
+ max_depth : int
+ The depth of the smallest HEALPix cells contained in the MOC.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC
+ """
+ hp = HEALPix(nside=(1 << max_depth), order='nested')
+ ipix = hp.lonlat_to_healpix(lon, lat)
+
+ shift = 2 * (MOCSkyRegion.HPY_MAX_DEPTH - max_depth)
+ itvs = np.vstack((ipix << shift, (ipix + 1) << shift)).T
+
+ itv = IntervalSet(itvs)
+ return cls(itv, meta, visual)
+
+ @classmethod
+ def from_polygon_skycoord(cls, skycoord, inside=None, max_depth=10, meta=None, visual=None):
+ """
+ Creates a MOC from a polygon.
+
+ The polygon is given as an `astropy.coordinates.SkyCoord` that contains the
+ vertices of the polygon. Concave and convex polygons are accepted but
+ self-intersecting ones are currently not properly handled.
+
+ Parameters
+ ----------
+ skycoord : `astropy.coordinates.SkyCoord`
+ The sky coordinates defining the vertices of a polygon. It can describe a convex or
+ concave polygon but not a self-intersecting one.
+ inside : `astropy.coordinates.SkyCoord`, optional
+ A point that will be inside the MOC is needed as it is not possible to determine the inside area of a polygon
+ on the unit sphere (there is no infinite area that can be considered as the outside because on the sphere,
+ a closed polygon delimits two finite areas).
+ Possible improvement: take the inside area as the one covering the smallest region on the sphere.
+
+ If inside=None (default behavior), the mean of all the vertices is taken as lying inside the polygon. That approach may not work for
+ concave polygons.
+ max_depth : int, optional
+ The resolution of the MOC. Set to 10 by default.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC
+ """
+ return MOCSkyRegion.from_polygon(lon=skycoord.icrs.ra, lat=skycoord.icrs.dec,
+ inside=inside, max_depth=max_depth, meta=meta, visual=visual)
+
+ @classmethod
+ def from_polygon(cls, lon, lat, inside=None, max_depth=10, meta=None, visual=None):
+ """
+ Creates a MOC from a polygon
+
+ The polygon is given as lon and lat `astropy.units.Quantity` that define the
+ vertices of the polygon. Concave and convex polygons are accepted but
+ self-intersecting ones are currently not properly handled.
+
+ Parameters
+ ----------
+ lon : `astropy.units.Quantity`
+ The longitudes defining the polygon. Can describe convex and
+ concave polygons but not self-intersecting ones.
+ lat : `astropy.units.Quantity`
+ The latitudes defining the polygon. Can describe convex and concave
+ polygons but not self-intersecting ones.
+ inside : `astropy.coordinates.SkyCoord`, optional
+ A point that will be inside the MOC is needed as it is not possible to determine the inside area of a polygon
+ on the unit sphere (there is no infinite area that can be considered as the outside because on the sphere,
+ a closed polygon delimits two finite areas).
+ Possible improvement: take the inside area as the one covering the smallest region on the sphere.
+
+ If inside=None (default behavior), the mean of all the vertices is taken as lying inside the polygon. That approach may not work for
+ concave polygons.
+ max_depth : int, optional
+ The resolution of the MOC. Set to 10 by default.
+ meta : `~regions.RegionMeta` object, optional
+ A dictionary which stores the meta attributes of this region.
+ visual : `~regions.RegionVisual` object, optional
+ A dictionary which stores the visual meta attributes of this region.
+
+ Returns
+ -------
+ result : `regions.MOCSkyRegion`
+ The resulting MOC
+ """
+ from .polygon import PolygonComputer
+
+ polygon_computer = PolygonComputer(lon, lat, inside, max_depth)
+ # Create the moc from the python dictionary
+
+ moc = MOCSkyRegion.from_json(polygon_computer.ipix, meta=meta, visual=visual)
+ # We degrade it to the user-requested order
+ if polygon_computer.degrade_to_max_depth:
+ moc = moc.degrade_to_depth(max_depth)
+
+ return moc
+
+ @property
+ def sky_fraction(self):
+ """
+ Sky fraction covered by the MOC
+ """
+ pix = self._best_res_pixels()
+ num_pix = pix.size
+ return num_pix / float(3 << (2*(self.max_depth + 1)))
diff --git a/regions/shapes/moc/interval_set.py b/regions/shapes/moc/interval_set.py
new file mode 100644
index 00000000..0a5a642d
--- /dev/null
+++ b/regions/shapes/moc/interval_set.py
@@ -0,0 +1,310 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function, unicode_literals
+from .py23_compat import range, int
+
+import copy
+import numpy as np
+
+from astropy_healpix import uniq_to_level_ipix
+
+class IntervalSet:
+ """
+ Internal data structure for storing a MOC.
+
+ Please refer to the doc example page. See :ref:`moc`.
+
+ Parameters
+ ----------
+ data : `~numpy.ndarray`
+ A N rows by 2 columns numpy array, each row representing one interval.
+ make_consistent : bool, optional
+ True by default. Remove the overlapping intervals from the data.
+ """
+ HPY_MAX_DEPTH = 29
+
+ def __init__(self, data=None, make_consistent=True):
+ data = np.array([], dtype=np.int64) if data is None else data
+ self._data = data
+ # Merging of the overlapping intervals
+ if make_consistent:
+ self._merge_intervals()
+
+ def copy(self):
+ """
+ Deepcopy of self.
+
+ Returns
+ -------
+ itv_s : `IntervalSet`
+ A copy of self
+ """
+ return copy.deepcopy(self)
+
+ def __repr__(self):
+ # Simply print the numpy array storing the intervals
+ return "{0}".format(self._data)
+
+ def __eq__(self, other):
+ """
+ Test equality of two `IntervalSet` objects
+
+ Parameters
+ ----------
+ other : `IntervalSet`
+ The other `IntervalSet` object
+
+ Returns
+ -------
+ res : bool
+ True if all the intervals from the two `IntervalSet` are equal one by one.
+ """
+ return np.all(self._data == other._data)
+
+ @property
+ def min(self):
+ return self._data.min()
+
+ @property
+ def max(self):
+ return self._data.max()
+
+ def empty(self):
+ """
+ Return True if there is no interval stored.
+ """
+ return self._data.size == 0
+
+ def _merge_intervals(self):
+ """
+ Merge overlapping intervals.
+
+ This method is called only one time at the end of the constructor.
+ """
+ res = []
+ start = stop = None
+ itvs_l = self._data.tolist()
+
+ for itv in sorted(itvs_l):
+ if start is None:
+ start, stop = itv
+ continue
+
+ # gap between intervals
+ if itv[0] > stop:
+ res.append((start, stop))
+ start, stop = itv
+ else:
+ # merge intervals
+ if itv[1] > stop:
+ stop = itv[1]
+
+ if start is not None and stop is not None:
+ res.append((start, stop))
+
+ self._data = np.asarray(res, dtype=np.int64)
+
+ def union(self, other):
+ """
+ Returns the union between two `IntervalSet` instances.
+
+ Parameters
+ ----------
+ other : `IntervalSet`
+ The other `IntervalSet` instance.
+ Returns
+ -------
+ res : `IntervalSet`
+ A new `IntervalSet` instance.
+ """
+ res = IntervalSet()
+ if other.empty():
+ res._data = self._data
+ elif self.empty():
+ res._data = other._data
+ else:
+ # res has no overlapping intervals
+ res._data = IntervalSet.merge(self._data,
+ other._data,
+ lambda in_a, in_b: in_a or in_b)
+ return res
+
+ def difference(self, other):
+ """
+ Returns the difference between two `IntervalSet` instances.
+
+ Parameters
+ ----------
+ other : `IntervalSet`
+ The other `IntervalSet` instance.
+ Returns
+ -------
+ res : `IntervalSet`
+ A new `IntervalSet` instance.
+ """
+ res = IntervalSet()
+ if other.empty() or self.empty():
+ res._data = self._data
+ else:
+ res._data = IntervalSet.merge(self._data,
+ other._data,
+ lambda in_a, in_b: in_a and not in_b)
+ return res
+
+ def intersection(self, other):
+ """
+ Return the intersection between self and ``another_is``.
+
+ Parameters
+ ----------
+ other : `IntervalSet`
+ The other `IntervalSet` object.
+ Returns
+ -------
+ res : `IntervalSet`
+ A new `IntervalSet` instance.
+ """
+ res = IntervalSet()
+ if not other.empty() and not self.empty():
+ res._data = IntervalSet.merge(self._data,
+ other._data,
+ lambda in_a, in_b: in_a and in_b)
+ return res
+
+ @classmethod
+ def to_uniq_itv_s(cls, itv_s):
+ """
+ Convert invervals (i.e. using the nested numbering scheme) to intervals of uniq numbers.
+
+ Parameters
+ ----------
+ itv_s : `IntervalSet`
+ A `IntervalSet` instance.
+
+ Returns
+ -------
+ result : `IntervalSet`
+ A new `IntervalSet` instance containing intervals in the uniq format.
+ """
+ r2 = itv_s.copy()
+ uniq_itvs_l = []
+
+ if r2.empty():
+ return IntervalSet()
+
+ depth = 0
+ while not r2.empty():
+ shift = int(2 * (IntervalSet.HPY_MAX_DEPTH - depth))
+ ofs = (int(1) << shift) - 1
+ ofs2 = int(1) << (2 * depth + 2)
+
+ r4 = []
+ for iv in r2._data:
+ a = (int(iv[0]) + ofs) >> shift
+ b = int(iv[1]) >> shift
+
+ c = a << shift
+ d = b << shift
+ if d > c:
+ r4.append((c, d))
+ uniq_itvs_l.append((a + ofs2, b + ofs2))
+
+ if len(r4) > 0:
+ r4_is = IntervalSet(np.asarray(r4, dtype=np.int64))
+ r2 = r2.difference(r4_is)
+
+ depth += 1
+
+ uniq_itvs = np.asarray(uniq_itvs_l, dtype=np.int64)
+ return IntervalSet(uniq_itvs)
+
+ @classmethod
+ def from_uniq_itv_s(cls, uniq_itv_s):
+ """
+ Convert uniq numbers intervals to nested intervals.
+
+ Parameters
+ ----------
+ uniq_itv_s : `IntervalSet`
+ A uniq formatted `IntervalSet` instance.
+
+ Returns
+ -------
+ itv_s : `IntervalSet`
+ A `IntervalSet` instance containing nested indices.
+ """
+ itv_s = IntervalSet()
+ # Appending a list is faster than appending a numpy array
+ rtmp = []
+ last_level = 0
+ uniq_itvs = uniq_itv_s._data
+ shift = 2 * IntervalSet.HPY_MAX_DEPTH
+ for uniq_itv in uniq_itvs:
+ for j in np.arange(uniq_itv[0], uniq_itv[1], dtype=np.int64):
+ level, ipix = uniq_to_level_ipix(j)
+
+ if level != last_level:
+ itv_s = itv_s.union(IntervalSet(np.asarray(rtmp, dtype=np.int64)))
+ rtmp = []
+ last_level = level
+ shift = 2 * (IntervalSet.HPY_MAX_DEPTH - level)
+
+ rtmp.append((ipix << shift, (ipix + 1) << shift))
+
+ itv_s = itv_s.union(IntervalSet(np.asarray(rtmp, dtype=np.int64)))
+ return itv_s
+
+ @staticmethod
+ def merge(a_itvs, b_itvs, op):
+ """
+ Merge two lists of intervals according to the boolean function op
+
+ ``a_itvs`` and ``b_itvs`` need to be sorted and consistent (i.e. no overlapping intervals).
+ This operation keeps the resulting `IntervalSet` sorted and consistent.
+
+ Parameters
+ ----------
+ a_itvs : `~numpy.ndarray`
+ A N rows by 2 columns numpy array storing non overlapping and sorted intervals.
+ b_itvs : `~numpy.ndarray`
+ A N rows by 2 columns numpy array storing non overlapping and sorted intervals.
+ op : `function`
+ Lambda function taking two params and returning the result of the operation between
+ these two params.
+ Exemple : lambda a, b: a and b describes the intersection of ``a_itvs`` and
+ ``b_itvs``.
+
+ Returns
+ -------
+ res : `~numpy.ndarray`
+ A N rows by 2 columns numpy array containing intervals resulting from the op between ``a_itvs`` and ``b_itvs``.
+ """
+ a_endpoints = a_itvs.flatten().tolist()
+ b_endpoints = b_itvs.flatten().tolist()
+
+ sentinel = max(a_endpoints[-1], b_endpoints[-1]) + 1
+
+ a_endpoints += [sentinel]
+ b_endpoints += [sentinel]
+
+ a_index = 0
+ b_index = 0
+
+ res = []
+
+ scan = min(a_endpoints[0], b_endpoints[0])
+ while scan < sentinel:
+ in_a = not ((scan < a_endpoints[a_index]) ^ (a_index % 2))
+ in_b = not ((scan < b_endpoints[b_index]) ^ (b_index % 2))
+ in_res = op(in_a, in_b)
+
+ if in_res ^ (len(res) % 2):
+ res += [scan]
+ if scan == a_endpoints[a_index]:
+ a_index += 1
+ if scan == b_endpoints[b_index]:
+ b_index += 1
+
+ scan = min(a_endpoints[a_index], b_endpoints[b_index])
+
+ # Reshape the resulting list of intervals to a Nx2 numpy array
+ return np.asarray(res, dtype=np.int64).reshape((-1, 2))
diff --git a/regions/shapes/moc/plot/__init__.py b/regions/shapes/moc/plot/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/regions/shapes/moc/plot/axis_viewport.py b/regions/shapes/moc/plot/axis_viewport.py
new file mode 100644
index 00000000..f799a877
--- /dev/null
+++ b/regions/shapes/moc/plot/axis_viewport.py
@@ -0,0 +1,11 @@
+def set(ax, wcs):
+ fig = ax.get_figure()
+ width_px, height_px = fig.get_size_inches() * fig.dpi
+
+ x_min = 0
+ x_max = width_px - 1
+ y_min = 0
+ y_max = height_px - 1
+
+ ax.set_xlim([x_min, x_max])
+ ax.set_ylim([y_min, y_max])
diff --git a/regions/shapes/moc/plot/border.py b/regions/shapes/moc/plot/border.py
new file mode 100644
index 00000000..44bb0cd7
--- /dev/null
+++ b/regions/shapes/moc/plot/border.py
@@ -0,0 +1,98 @@
+import numpy as np
+
+from astropy_healpix import HEALPix
+from astropy.coordinates import ICRS
+
+from astropy.wcs.utils import skycoord_to_pixel
+
+from astropy_healpix import level_to_nside, nside_to_npix
+
+def border(moc, ax, wcs, **kw_mpl_pathpatch):
+ from matplotlib.path import Path
+ from matplotlib.patches import PathPatch
+
+ from .utils import build_plotting_moc
+ moc_to_plot = build_plotting_moc(moc, wcs)
+
+ if moc_to_plot.empty():
+ return
+
+ max_depth = moc_to_plot.max_depth
+ nside = level_to_nside(max_depth)
+ hp = HEALPix(nside=nside, order='nested', frame=ICRS())
+ ipixels_open = moc_to_plot._best_res_pixels()
+
+ # Take the complement if the MOC covers more than half of the sky
+ num_ipixels = nside_to_npix(nside)
+ sky_fraction = ipixels_open.shape[0] / float(num_ipixels)
+
+ if sky_fraction > 0.5:
+ ipixels_all = np.arange(num_ipixels, dtype=np.int64)
+ ipixels_open = np.setdiff1d(ipixels_all, ipixels_open, assume_unique=True)
+
+ ipix_neigh = hp.neighbours(ipixels_open)
+ # Select the direct neighbors (i.e. those in WEST, NORTH, EAST and SOUTH directions)
+ ipix_neigh = ipix_neigh[[0, 2, 4, 6], :]
+
+ r1 = np.in1d(ipix_neigh[0, :], ipixels_open)
+ r2 = np.in1d(ipix_neigh[1, :], ipixels_open)
+ r3 = np.in1d(ipix_neigh[2, :], ipixels_open)
+ r4 = np.in1d(ipix_neigh[3, :], ipixels_open)
+
+ ipix_moc = np.vstack((r1, r2, r3, r4))
+
+ west_edge = ipix_moc[0, :]
+ south_edge = ipix_moc[1, :]
+ east_edge = ipix_moc[2, :]
+ north_edge = ipix_moc[3, :]
+
+ num_ipix_moc = ipix_moc.sum(axis=0)
+
+ ipixels_border_id = (num_ipix_moc < 4)
+ # The border of each HEALPix cells is drawn one at a time
+ path_vertices_l = []
+ codes = []
+
+ west_border = west_edge[ipixels_border_id]
+ south_border = south_edge[ipixels_border_id]
+ east_border = east_edge[ipixels_border_id]
+ north_border = north_edge[ipixels_border_id]
+ ipixels_border = ipixels_open[ipixels_border_id]
+ ipix_boundaries = hp.boundaries_skycoord(ipixels_border, step=1)
+
+ # Projection on the given WCS
+ xp, yp = skycoord_to_pixel(coords=ipix_boundaries, wcs=wcs)
+ from . import culling_backfacing_cells
+ xp, yp, frontface_id = culling_backfacing_cells.backface_culling(xp, yp)
+
+ west_border = west_border[frontface_id]
+ south_border = south_border[frontface_id]
+ east_border = east_border[frontface_id]
+ north_border = north_border[frontface_id]
+
+ for i in range(xp.shape[0]):
+ vx = xp[i]
+ vy = yp[i]
+ if not south_border[i]:
+ path_vertices_l += [(vx[0], vy[0]), (vx[1], vy[1]), (0, 0)]
+ codes += [Path.MOVETO] + [Path.LINETO] + [Path.CLOSEPOLY]
+
+ if not west_border[i]:
+ path_vertices_l += [(vx[1], vy[1]), (vx[2], vy[2]), (0, 0)]
+ codes += [Path.MOVETO] + [Path.LINETO] + [Path.CLOSEPOLY]
+
+ if not north_border[i]:
+ path_vertices_l += [(vx[2], vy[2]), (vx[3], vy[3]), (0, 0)]
+ codes += [Path.MOVETO] + [Path.LINETO] + [Path.CLOSEPOLY]
+
+ if not east_border[i]:
+ path_vertices_l += [(vx[3], vy[3]), (vx[0], vy[0]), (0, 0)]
+ codes += [Path.MOVETO] + [Path.LINETO] + [Path.CLOSEPOLY]
+
+ path = Path(path_vertices_l, codes)
+ perimeter_patch = PathPatch(path, **kw_mpl_pathpatch)
+
+ ax.add_patch(perimeter_patch)
+
+ from . import axis_viewport
+ axis_viewport.set(ax, wcs)
diff --git a/regions/shapes/moc/plot/culling_backfacing_cells.py b/regions/shapes/moc/plot/culling_backfacing_cells.py
new file mode 100644
index 00000000..c23d4811
--- /dev/null
+++ b/regions/shapes/moc/plot/culling_backfacing_cells.py
@@ -0,0 +1,94 @@
+import numpy as np
+
+from astropy_healpix import HEALPix
+from astropy.coordinates import ICRS
+from astropy.wcs.utils import skycoord_to_pixel
+
+def backface_culling(xp, yp):
+ # Remove cells crossing the MOC after projection
+ # The remaining HEALPix cells are used for computing the patch of the MOC
+ vx = xp
+ vy = yp
+
+ def cross_product(vx, vy, i):
+ cur = i
+ prev = (i - 1) % 4
+ next = (i + 1) % 4
+
+ # Construct the first vector from A to B
+ x1 = vx[:, cur] - vx[:, prev]
+ y1 = vy[:, cur] - vy[:, prev]
+ z1 = np.zeros(x1.shape)
+
+ v1 = np.vstack((x1, y1, z1)).T
+ # Construct the second vector from B to C
+ x2 = vx[:, next] - vx[:, cur]
+ y2 = vy[:, next] - vy[:, cur]
+ z2 = np.zeros(x2.shape)
+
+ v2 = np.vstack((x2, y2, z2)).T
+ # Compute the cross product between the two
+ return np.cross(v1, v2)
+
+ # A ----- B
+ # \ |
+ # D-----C
+ # Compute the cross product between AB and BC
+ # and the cross product between BC and CD
+ ABC = cross_product(vx, vy, 1)
+ CDA = cross_product(vx, vy, 3)
+
+ frontface_cells = (ABC[:, 2] < 0) & (CDA[:, 2] < 0)
+ frontface_cells = np.asarray(frontface_cells)
+
+ vx = vx[frontface_cells]
+ vy = vy[frontface_cells]
+
+ return vx, vy, frontface_cells
+
+def from_moc(depth_ipix_d, wcs):
+ # Create a new MOC that do not contain the HEALPix
+ # cells that are backfacing the projection
+ depths = [int(depth_str) for depth_str in depth_ipix_d.keys()]
+ min_depth = min(depths)
+ max_depth = max(depths)
+ ipixels = np.asarray(depth_ipix_d[str(min_depth)])
+
+ max_split_depth = max_depth
+
+ ipix_d = {}
+ for depth in range(min_depth, max_split_depth + 1):
+ hp = HEALPix(nside=(1 << depth), order='nested', frame=ICRS())
+
+ ipix_boundaries = hp.boundaries_skycoord(ipixels, step=1)
+ # Projection on the given WCS
+ xp, yp = skycoord_to_pixel(coords=ipix_boundaries, wcs=wcs)
+ xp, yp, frontface_id = backface_culling(xp, yp)
+
+ # Get the pixels which are backfacing the projection
+ backfacing_ipix = ipixels[~frontface_id]
+ frontface_ipix = ipixels[frontface_id]
+
+ depth_str = str(depth)
+ ipix_d.update({depth_str: frontface_ipix})
+
+ if depth < max_depth:
+ next_depth = str(depth + 1)
+ ipixels = []
+ if next_depth in depth_ipix_d:
+ ipixels = depth_ipix_d[next_depth]
+
+ for bf_ipix in backfacing_ipix:
+ child_bf_ipix = bf_ipix << 2
+ ipixels.extend([child_bf_ipix,
+ child_bf_ipix + 1,
+ child_bf_ipix + 2,
+ child_bf_ipix + 3])
+
+ ipixels = np.asarray(ipixels)
+
+ for depth in range(max_split_depth+1, max_depth+1):
+ depth_str = str(depth)
+ ipix_d.update({depth_str: depth_ipix_d[depth_str]})
+
+ return ipix_d
diff --git a/regions/shapes/moc/plot/fill.py b/regions/shapes/moc/plot/fill.py
new file mode 100644
index 00000000..78f1b231
--- /dev/null
+++ b/regions/shapes/moc/plot/fill.py
@@ -0,0 +1,106 @@
+import numpy as np
+
+from astropy_healpix import HEALPix
+from astropy.coordinates import ICRS
+
+from astropy.wcs.utils import skycoord_to_pixel
+
+from .utils import build_plotting_moc
+from . import culling_backfacing_cells
+
+from astropy_healpix import level_to_nside
+
+def compute_healpix_vertices(depth, ipix, wcs, origin):
+ from matplotlib.path import Path
+
+ path_vertices = np.array([])
+ codes = np.array([])
+
+ depth = int(depth)
+
+ step = 1
+ if depth < 3:
+ step = 2
+
+ nside = level_to_nside(depth)
+ hp = HEALPix(order="nested", nside=nside, frame=ICRS())
+ ipix_boundaries = hp.boundaries_skycoord(ipix, step=step)
+ # Projection on the given WCS
+ xp, yp = skycoord_to_pixel(ipix_boundaries, wcs=wcs)
+
+ # Add offset origin given by the user (default is origin=(0, 0))
+ xp += origin[0]
+ yp += origin[1]
+
+ c1 = np.vstack((xp[:, 0], yp[:, 0])).T
+ c2 = np.vstack((xp[:, 1], yp[:, 1])).T
+ c3 = np.vstack((xp[:, 2], yp[:, 2])).T
+ c4 = np.vstack((xp[:, 3], yp[:, 3])).T
+
+ if depth < 3:
+ c5 = np.vstack((xp[:, 4], yp[:, 4])).T
+ c6 = np.vstack((xp[:, 5], yp[:, 5])).T
+ c7 = np.vstack((xp[:, 6], yp[:, 6])).T
+ c8 = np.vstack((xp[:, 7], yp[:, 7])).T
+
+ cells = np.hstack((c1, c2, c3, c4, c5, c6, c7, c8, np.zeros((c1.shape[0], 2))))
+
+ path_vertices = cells.reshape((9*c1.shape[0], 2))
+ single_code = np.array([Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY])
+ else:
+ cells = np.hstack((c1, c2, c3, c4, np.zeros((c1.shape[0], 2))))
+
+ path_vertices = cells.reshape((5*c1.shape[0], 2))
+ single_code = np.array([Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY])
+
+ codes = np.tile(single_code, c1.shape[0])
+
+ return path_vertices, codes
+
+def compute_the_patches(moc, wcs, origin=(0, 0)):
+ depth_ipix_d = moc.serialize(format="json")
+ depth_ipix_clean_d = culling_backfacing_cells.from_moc(depth_ipix_d=depth_ipix_d, wcs=wcs)
+
+ patches = []
+ for depth, ipix in depth_ipix_clean_d.items():
+ patch = compute_healpix_vertices(depth=depth,
+ ipix=ipix,
+ wcs=wcs,
+ origin=origin)
+ patches.append(patch)
+
+ return depth_ipix_clean_d.keys(), patches
+
+def build_mpl_pathpatch(patches, wcs, **kw_mpl_pathpatch):
+ from matplotlib.path import Path
+ from matplotlib.patches import PathPatch
+
+ first_patch = patches[0]
+ vertices_first_patch, codes_first_patch = first_patch
+ path_vertices = np.array(vertices_first_patch)
+ path_codes = np.array(codes_first_patch)
+
+ for vertices, codes in patches[1:]:
+ path_vertices = np.vstack((path_vertices, vertices))
+ path_codes = np.hstack((path_codes, codes))
+
+ path = Path(path_vertices, path_codes)
+ pathpatch_mpl = PathPatch(path, **kw_mpl_pathpatch)
+
+ return pathpatch_mpl
+
+def fill(moc, wcs, origin, **kw_mpl_pathpatch):
+ # Simplify the MOC for plotting purposes:
+ # 1. Degrade the MOC if the FOV is enough big so that we cannot see the smallest HEALPix cells.
+ # 2. For small FOVs, plot the MOC & POLYGONAL_MOC_FROM_FOV.
+ moc_to_plot = build_plotting_moc(moc=moc, wcs=wcs)
+
+ # If the FOV contains no cells, then moc_to_plot (i.e. the intersection between the moc
+ # and the MOC created from the FOV polygon) will be empty.
+ # If it is the case, we exit the method without doing anything.
+ if moc_to_plot.empty():
+ return None
+
+ _, patches = compute_the_patches(moc=moc_to_plot, wcs=wcs, origin=origin)
+ pathpatch_mpl = build_mpl_pathpatch(patches=patches, wcs=wcs, **kw_mpl_pathpatch)
+ return pathpatch_mpl
diff --git a/regions/shapes/moc/plot/utils.py b/regions/shapes/moc/plot/utils.py
new file mode 100644
index 00000000..099671c7
--- /dev/null
+++ b/regions/shapes/moc/plot/utils.py
@@ -0,0 +1,64 @@
+import numpy as np
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+from astropy.wcs.utils import pixel_to_skycoord
+
+import warnings
+
+def build_plotting_moc(moc, wcs):
+ # Get the WCS cdelt giving the deg.px^(-1) resolution.
+ cdelt = wcs.wcs.cdelt
+ # Convert in rad.px^(-1)
+ cdelt = np.abs((2*np.pi/360)*cdelt[0])
+ # Get the minimum depth such as the resolution of a cell is contained in 1px.
+ depth_res = int(np.floor(np.log2(np.sqrt(np.pi/3)/cdelt)))
+ depth_res = max(depth_res, 0)
+ # Degrade the moc to that depth for plotting purposes. It is not necessary to plot pixels
+ # that we will not see because they are contained in 1px.
+ moc_plot = moc
+ if moc.max_depth > depth_res:
+ moc_plot = moc.degrade_to_depth(depth_res)
+
+ # Get the MOC delimiting the FOV polygon
+ width_px = int(wcs.wcs.crpix[0]*2.) # Supposing the wcs is centered in the axis
+ heigth_px = int(wcs.wcs.crpix[1]*2.)
+
+ # Compute the sky coordinate path delimiting the viewport.
+ # It consists of a closed polygon of (4 - 1)*4 = 12 vertices
+ x_px = np.linspace(0, width_px, 4)
+ y_px = np.linspace(0, heigth_px, 4)
+
+ X, Y = np.meshgrid(x_px, y_px)
+
+ X_px = np.append(X[0, :-1], X[:-1, -1])
+ X_px = np.append(X_px, X[-1, 1:][::-1])
+ X_px = np.append(X_px, X[:-1, 0])
+
+ Y_px = np.append(Y[0, :-1], Y[:-1, -1])
+ Y_px = np.append(Y_px, Y[-1, :-1])
+ Y_px = np.append(Y_px, Y[1:, 0][::-1])
+
+ # Disable the output of warnings when encoutering NaNs.
+ warnings.filterwarnings("ignore")
+ # Inverse projection from pixel coordinate space to the world coordinate space
+ viewport = pixel_to_skycoord(X_px, Y_px, wcs)
+ # If one coordinate is a NaN we exit the function and do not go further
+ ra_deg, dec_deg = viewport.icrs.ra.deg, viewport.icrs.dec.deg
+ warnings.filterwarnings("default")
+
+ if np.isnan(ra_deg).any() or np.isnan(dec_deg).any():
+ return moc_plot
+
+ center_x_px, center_y_px = wcs.wcs.crpix[0], wcs.wcs.crpix[1]
+ inside = pixel_to_skycoord(center_x_px, center_y_px, wcs)
+
+ # Import MOC here to avoid circular imports
+ from ..core import MOCSkyRegion
+ # Create a rough MOC (depth=3 is sufficient) from the viewport
+ moc_viewport = MOCSkyRegion.from_polygon_skycoord(viewport, max_depth=3, inside=inside)
+
+ # The moc to plot is the INPUT_MOC & MOC_VIEWPORT. For small FOVs this can reduce
+ # a lot the time to draw the MOC along with its borders.
+ moc_plot = moc_plot.intersection(moc_viewport)
+ return moc_plot
diff --git a/regions/shapes/moc/plot/wcs.py b/regions/shapes/moc/plot/wcs.py
new file mode 100644
index 00000000..dabe8d18
--- /dev/null
+++ b/regions/shapes/moc/plot/wcs.py
@@ -0,0 +1,94 @@
+import numpy as np
+
+from astropy import coordinates
+from astropy import wcs
+
+import astropy.units as u
+
+class WCS:
+ """
+ Creates a WCS for vizualizing a `~regions.MOCSkyRegion` in a matplotlib axis.
+
+ Parameters
+ ----------
+ fig : `~matplotlib.pyplot.figure`
+ The matplotlib figure used for plotting the MOC.
+ fov : `~astropy.units.Quantity`
+ Size of the field of view.
+ center : `~astropy.coordinates.SkyCoord`, optional
+ World coordinates matching with the center of the plot. Default to (0 deg, 0 deg) (in ICRS frame).
+ coordsys : str, optional
+ Coordinate system. Default to "icrs". Must be in ["icrs", "galactic"].
+ projection : str, optional
+ World base -> Image base projection type. See http://docs.astropy.org/en/stable/wcs/#supported-projections for
+ the projections currently supported in astropy. Default to Aitoff.
+ rotation : `~astropy.coordinates.Angle`, optional
+ The angle of rotation. Default to no rotation.
+
+ Returns
+ -------
+ wcs : `~astropy.wcs.WCS`
+ The WCS that can be passed to `regions.MOCSkyRegion.fill` and `regions.MOCSkyRegion.border`.
+
+ Examples
+ --------
+ >>> 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
+ >>> 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=(15, 15))
+ >>> # Define a WCS as a context
+ >>> with WCS(fig,
+ ... fov=200 * u.deg,
+ ... center=SkyCoord(0, 20, 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.
+ ... 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.grid(color="black", linestyle="dotted")
+ """
+ def __init__(self,
+ fig,
+ fov,
+ center=coordinates.SkyCoord(0, 0, unit="deg", frame="icrs"),
+ coordsys="icrs",
+ projection="AIT",
+ rotation=coordinates.Angle(0, u.radian)):
+ self.w = wcs.WCS(naxis=2)
+
+ width_px, height_px = fig.get_size_inches() * float(fig.dpi)
+
+ cdelt_x = fov.to_value("deg")/float(width_px)
+ cdelt_y = fov.to_value("deg")/float(height_px)
+
+ self.w.wcs.crpix = [width_px/2.0, height_px/2.0]
+ self.w.wcs.cdelt = [-cdelt_x, cdelt_x]
+
+ if coordsys == 'icrs':
+ self.w.wcs.crval = [center.icrs.ra.deg, center.icrs.dec.deg]
+ self.w.wcs.ctype = ['RA---' + projection, 'DEC--' + projection]
+ elif coordsys == 'galactic':
+ self.w.wcs.crval = [center.galactic.l.deg, center.galactic.b.deg]
+ self.w.wcs.ctype = ['GLON-' + projection, 'GLAT-' + projection]
+
+ theta = rotation.radian
+ self.w.wcs.pc = [
+ [np.cos(theta), -np.sin(theta)],
+ [np.sin(theta), np.cos(theta)],
+ ]
+
+ def __enter__(self):
+ return self.w
+
+ def __exit__(self, exception_type, exception_value, traceback):
+ pass
diff --git a/regions/shapes/moc/polygon.py b/regions/shapes/moc/polygon.py
new file mode 100644
index 00000000..69ad4dce
--- /dev/null
+++ b/regions/shapes/moc/polygon.py
@@ -0,0 +1,217 @@
+import numpy as np
+
+from astropy import units as u
+from astropy.coordinates import ICRS, Galactic
+from astropy.coordinates import SkyCoord
+
+from astropy_healpix import HEALPix
+from astropy_healpix.core import boundaries_lonlat
+
+try:
+ from spherical_geometry.polygon import SphericalPolygon
+ from spherical_geometry import great_circle_arc
+ from spherical_geometry import vector
+except ImportError as e:
+ print("Unable to find spherical-geometry package installed. See https://github.com/spacetelescope/spherical_geometry")
+ raise
+
+class PolygonComputer:
+ def polygon_contains_ipix(self, ipix):
+ def poly_contains_ipix_vertices(poly, ipix):
+ # Get the ipix polygon vertices as a Nx3 numpy array
+ # Remove the last vertex as it counts double (closed polygon)
+ ipix_vertices = list(ipix.points)[0][:-1, :]
+
+ for ipix_vertex in ipix_vertices:
+ if not poly.contains_point(ipix_vertex):
+ return False
+
+ return True
+
+ def poly_crossing_ipix(poly, ipix):
+ # ipix and poly are spherical shapes
+ poly_points = list(poly.points)[0]
+ x_poly, y_poly, z_poly = (poly_points[:, 0], poly_points[:, 1], poly_points[:, 2])
+
+ ipix_points = list(ipix.points)[0]
+ x_ipix, y_ipix, z_ipix = (ipix_points[:, 0], ipix_points[:, 1], ipix_points[:, 2])
+
+ A_poly = np.stack((x_poly, y_poly, z_poly)).T
+ B_poly = A_poly
+ B_poly = np.append(B_poly, [B_poly[1]], axis=0)
+ B_poly = B_poly[1:]
+
+ A_poly = A_poly[:-1]
+ B_poly = B_poly[:-1]
+
+ for i in range(len(ipix_points) - 1):
+ A_ipix = (x_ipix[i], y_ipix[i], z_ipix[i])
+ B_ipix = (x_ipix[i+1], y_ipix[i+1], z_ipix[i+1])
+
+ inter = great_circle_arc.intersects(A_poly, B_poly, A_ipix, B_ipix)
+ if inter.any():
+ return True
+ return False
+
+ return poly_contains_ipix_vertices(self.polygon, ipix) and \
+ (not poly_contains_ipix_vertices(ipix, self.polygon)) and \
+ (not poly_crossing_ipix(self.polygon, ipix))
+
+ def _get_starting_depth(self):
+ def compute_angular_distance(n1, n2):
+ return np.arctan(np.linalg.norm(np.cross(n1, n2))/np.dot(n1, n2))
+
+ def to_xyz(lon, lat):
+ x = np.cos(lon) * np.cos(lat)
+ y = np.cos(lat) * np.sin(lon)
+ z = np.sin(lat)
+
+ return np.array([x, y, z], dtype=np.float64)
+
+ def max_distance_center_to_vertex(depth):
+ nside = 1 << depth
+
+ lat1 = np.arcsin(2 / 3.0)
+ lat2 = np.arcsin(1 - ((1 - 1.0/nside)**2 / 3.0))
+ lon1 = np.pi/(4 * nside)
+ lon2 = 0
+
+ # Convert lon, lat to xyz, vector
+ n1 = to_xyz(lon=lon1, lat=lat1)
+ n2 = to_xyz(lon=lon2, lat=lat2)
+
+ dist = compute_angular_distance(n1, n2)
+ return dist # in rad
+
+ # Get the polygon vertices as a Nx3 numpy array
+ # Remove the last vertex as it counts double (closed polygon)
+ p_vertices = np.asarray(list(self.polygon.points))[0][:-1, :]
+ # Get the center formed by the vertices
+ center = p_vertices.mean(axis=0)
+
+ # Normalize it so that it lies on the unit sphere
+ vector.normalize_vector(center, output=center)
+ center = np.asarray(center)
+ # Compute the maximum angular distance between the polygon vertices
+ # and its center. This refers to the Vector version of the
+ # Great-circle distance computation.
+ #
+ # See https://en.wikipedia.org/wiki/Great-circle_distance
+ max_d = -1
+
+ # Check if the polygon covers more than an hemisphere
+ covers_more_than_one_hemisphere = (self.polygon.area() > 2 * np.pi)
+
+ for vertex in p_vertices:
+ d = compute_angular_distance(center, vertex)
+ if covers_more_than_one_hemisphere:
+ d = np.pi - d
+
+ if d > max_d:
+ max_d = d
+
+ # Return the min depth so that max_d > max_center_to_vertex_ipix(depth)
+ depth = 0
+ while max_distance_center_to_vertex(depth) >= max_d:
+ depth = depth + 1
+
+ # Get the ipixels from astropy_healpix covering the cone of (center, radius) = (center, max_d)
+ lon_center, lat_center = vector.vector_to_lonlat(x=center[0], y=center[1], z=center[2], degrees=False)
+ hp = HEALPix(nside=(1 << depth), order='nested', frame=ICRS())
+ starting_iter_ipix = hp.cone_search_lonlat(lon=lon_center * u.rad, lat=lat_center * u.rad, radius=max_d * u.rad)
+
+ return depth, starting_iter_ipix
+
+ def _closes_numpy_2d_array(self, x):
+ tmp = np.zeros((x.shape[0], x.shape[1] + 1))
+ tmp[:, :-1] = x
+ tmp[:, -1] = x[:, 0]
+ return tmp
+
+ @property
+ def ipix(self):
+ return self.ipix_d
+
+ def __init__(self, ra, dec, inside=None, max_depth=10):
+ ra = ra.to(u.rad).value
+ dec = dec.to(u.rad).value
+ # Check if the vertices form a closed polygon
+ if ra[0] != ra[-1] or dec[0] != dec[-1]:
+ # If not, append the first vertex to ``vertices``
+ ra = np.append(ra, ra[0])
+ dec = np.append(dec, dec[0])
+ vertices = SkyCoord(ra=ra, dec=dec, unit="rad", frame="icrs")
+
+ if inside:
+ # Convert it to (x, y, z) cartesian coordinates on the sphere
+ inside = (inside.icrs.ra.rad, inside.icrs.dec.rad)
+
+ self.polygon = SphericalPolygon.from_lonlat(lon=ra, lat=dec, center=inside, degrees=False)
+
+ start_depth, ipixels = self._get_starting_depth()
+ end_depth = max_depth
+
+ # When the start depth returned is > to the depth requested
+ # For that specific case, we only do one iteration at start_depth
+ # Thus the MOC will contain the partially intersecting cells with the
+ # contained ones at start_depth
+
+ # And we degrade the MOC to the max_depth
+ self.degrade_to_max_depth = False
+ if start_depth > end_depth:
+ end_depth = start_depth
+ self.degrade_to_max_depth = True
+
+ self.ipix_d = {str(order): [] for order in range(start_depth, end_depth + 1)}
+
+ ## Iterative version of the algorithm: seems a bit faster than the recursive one
+ for depth in range(start_depth, end_depth + 1):
+ # Define a HEALPix at the current depth
+ hp = HEALPix(nside=(1 << depth), order='nested', frame=ICRS())
+
+ # Get the lon and lat of the corners of the pixels
+ # intersecting the polygon
+ lon, lat = hp.boundaries_lonlat(ipixels, step=1)
+ lon = lon.to(u.rad).value
+ lat = lat.to(u.rad).value
+
+ # closes the lon and lat array so that their first and last value matches
+ lon = self._closes_numpy_2d_array(lon)
+ lat = self._closes_numpy_2d_array(lat)
+
+ num_ipix_inter_poly = ipixels.shape[0]
+
+ # Define a 3d numpy array containing the corners coordinates of the intersecting pixels
+ # The first dim is the num of ipixels
+ # The second is the number of coordinates (5 as it defines the closed polygon of a HEALPix cell)
+ # The last is of size 2 (lon and lat)
+ shapes = np.vstack((lon.ravel(), lat.ravel())).T.reshape(num_ipix_inter_poly, 5, -1)
+
+ ipix_in_polygon_l = []
+ ipix_inter_polygon_l = []
+
+ for i in range(num_ipix_inter_poly):
+ shape = shapes[i]
+ # Definition of a SphericalPolygon from the border coordinates of a HEALPix cell
+ ipix_shape = SphericalPolygon.from_radec(lon=shape[:, 0], lat=shape[:, 1], degrees=False)
+ ipix = ipixels[i]
+
+ if self.polygon.intersects_poly(ipix_shape):
+ # If we are at the max depth then we direcly add to the MOC the intersecting ipixels
+ if depth == end_depth:
+ ipix_in_polygon_l.append(ipix)
+ else:
+ # Check whether polygon contains ipix or not
+ if self.polygon_contains_ipix(ipix_shape):
+ ipix_in_polygon_l.append(ipix)
+ else:
+ # The ipix is just intersecting without being contained in the polygon
+ # We split it in its 4 children
+ child_ipix = ipix << 2
+ ipix_inter_polygon_l.extend([child_ipix,
+ child_ipix + 1,
+ child_ipix + 2,
+ child_ipix + 3])
+
+ self.ipix_d.update({str(depth): ipix_in_polygon_l})
+ ipixels = np.asarray(ipix_inter_polygon_l, dtype=np.int64)
diff --git a/regions/shapes/moc/py23_compat.py b/regions/shapes/moc/py23_compat.py
new file mode 100644
index 00000000..d99e199b
--- /dev/null
+++ b/regions/shapes/moc/py23_compat.py
@@ -0,0 +1,19 @@
+"""
+py23_compat.py
+
+Python 2 / 3 compatibility layer.
+"""
+import sys
+
+__all__ = [
+ 'int', 'range',
+]
+
+PY2 = (sys.version_info.major == 2)
+
+if PY2:
+ int = long
+ range = xrange
+else:
+ int = int
+ range = range
diff --git a/regions/shapes/moc/utils.py b/regions/shapes/moc/utils.py
new file mode 100644
index 00000000..da495837
--- /dev/null
+++ b/regions/shapes/moc/utils.py
@@ -0,0 +1,26 @@
+def trailing_zeros(x):
+ bits = 0
+ # convention for x == 0 => return 0
+ if x == 0:
+ return 0
+
+ if not (x & 0xFFFFFFFF):
+ bits += 32
+ x >>= 32
+ if not (x & 0xFFFF):
+ bits += 16
+ x >>= 16
+ if not (x & 0xFF):
+ bits += 8
+ x >>= 8
+ if not (x & 0xF):
+ bits += 4
+ x >>= 4
+ if not (x & 0x3):
+ bits += 2
+ x >>= 2
+ if not (x & 1):
+ bits += 1
+ x >>= 1
+
+ return bits
diff --git a/regions/shapes/tests/data/P-GALEXGR6-AIS-FUV.fits b/regions/shapes/tests/data/P-GALEXGR6-AIS-FUV.fits
new file mode 100644
index 00000000..4418cae4
Binary files /dev/null and b/regions/shapes/tests/data/P-GALEXGR6-AIS-FUV.fits differ
diff --git a/regions/shapes/tests/data/P-SDSS9-r.fits b/regions/shapes/tests/data/P-SDSS9-r.fits
new file mode 100644
index 00000000..c9654df7
Binary files /dev/null and b/regions/shapes/tests/data/P-SDSS9-r.fits differ
diff --git a/regions/shapes/tests/test_interval_set.py b/regions/shapes/tests/test_interval_set.py
new file mode 100644
index 00000000..46bdd7aa
--- /dev/null
+++ b/regions/shapes/tests/test_interval_set.py
@@ -0,0 +1,85 @@
+import pytest
+import numpy as np
+from ..moc.interval_set import IntervalSet
+
+
+@pytest.fixture()
+def isets():
+ a = IntervalSet(np.asarray([(49, 73), (53, 54), (33, 63), (65, 80), (51, 80), (100, 126), (38, 68), (61, 72), (74, 102), (27, 43)]))
+ b = IntervalSet(np.asarray([(17, 26), (17, 41), (12, 31), (32, 61), (68, 90), (77, 105), (18, 27), (12, 35), (9, 37), (87, 97)]))
+ return dict(a=a, b=b)
+
+
+def test_interval_set_consistency(isets):
+ assert isets['a'] == IntervalSet(np.asarray([(27, 126)]))
+ assert isets['b'] == IntervalSet(np.asarray([(9, 61), (68, 105)]))
+
+
+def test_interval_set_union(isets):
+ assert isets['a'].union(isets['b']) == IntervalSet(np.array([(9, 126)]))
+ assert isets['a'].union(IntervalSet()) == IntervalSet(np.array([(27, 126)]))
+ assert IntervalSet().union(isets['a']) == IntervalSet(np.array([(27, 126)]))
+
+
+def test_interval_set_intersection(isets):
+ assert isets['a'].intersection(isets['b']) == IntervalSet(np.asarray([(27, 61), (68, 105)]))
+ assert isets['a'].intersection(IntervalSet()) == IntervalSet()
+ assert IntervalSet().intersection(isets['a']) == IntervalSet()
+
+
+def test_interval_set_difference(isets):
+ assert isets['a'].difference(isets['b']) == IntervalSet(np.asarray([(61, 68), (105, 126)]))
+ assert isets['b'].difference(isets['a']) == IntervalSet(np.asarray([(9, 27)]))
+ assert IntervalSet().difference(isets['a']) == IntervalSet()
+ assert isets['a'].difference(IntervalSet()) == isets['a']
+
+
+@pytest.fixture()
+def isets2():
+ nested1 = IntervalSet(np.array([[0, 1]]))
+ nuniq1 = IntervalSet(np.array([[4*4**29, 4*4**29 + 1]]))
+ nested2 = IntervalSet(np.array([[7, 76]]))
+ nuniq2 = IntervalSet(np.array([[1 + 4*4**27, 4 + 4*4**27],
+ [2 + 4*4**28, 4 + 4*4**28],
+ [16 + 4*4**28, 19 + 4*4**28],
+ [7 + 4*4**29, 8 + 4*4**29]]))
+ return dict(nested1=nested1, nuniq1=nuniq1,
+ nested2=nested2, nuniq2=nuniq2)
+
+
+def test_to_nuinq_interval_set(isets2):
+ assert IntervalSet.to_uniq_itv_s(isets2['nested1']) == isets2['nuniq1']
+ assert IntervalSet.to_uniq_itv_s(isets2['nested2']) == isets2['nuniq2']
+ # empty nested interval set
+ assert IntervalSet.to_uniq_itv_s(IntervalSet()) == IntervalSet()
+
+
+def test_from_nuinq_interval_set(isets2):
+ assert IntervalSet.from_uniq_itv_s(isets2['nuniq1']) == isets2['nested1']
+ assert IntervalSet.from_uniq_itv_s(isets2['nuniq2']) == isets2['nested2']
+ # empty nuniq interval set
+ assert IntervalSet.from_uniq_itv_s(IntervalSet()) == IntervalSet()
+
+
+def test_from_to_interval_set(isets2):
+ assert IntervalSet.from_uniq_itv_s(
+ IntervalSet.to_uniq_itv_s(isets2['nested1'])
+ ) == isets2['nested1']
+
+
+def test_interval_set_min(isets):
+ assert isets['a'].min == 27
+ assert isets['b'].min == 9
+ assert isets['a'].union(isets['b']).min == 9
+
+
+def test_interval_set_max(isets):
+ assert isets['a'].max == 126
+ assert isets['b'].max == 105
+ assert isets['a'].union(isets['b']).max == 126
+
+
+def test_repr_interval_set(isets):
+ assert repr(isets['a']) == "[[ 27 126]]"
+ assert repr(isets['b']) == "[[ 9 61]\n" \
+ " [ 68 105]]"
diff --git a/regions/shapes/tests/test_moc.py b/regions/shapes/tests/test_moc.py
new file mode 100644
index 00000000..deec416d
--- /dev/null
+++ b/regions/shapes/tests/test_moc.py
@@ -0,0 +1,268 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function
+
+import pytest
+import tempfile
+import numpy as np
+import copy
+import os
+
+from astropy.utils.data import get_pkg_data_filename
+from astropy.coordinates import SkyCoord, Angle
+import astropy.units as u
+from astropy.coordinates import ICRS
+from astropy.io import fits
+from astropy_healpix import HEALPix
+
+from ...core import PixCoord, RegionMask
+from ..._utils.examples import make_example_dataset
+from ...tests.helpers import make_simple_wcs
+
+from ..moc import MOCSkyRegion, WCS
+
+from .utils import HAS_MATPLOTLIB
+
+@pytest.fixture()
+def random_skycoords():
+ def generate(size):
+ return SkyCoord(np.random.uniform(0, 360, size),
+ np.random.uniform(-89, 89, size),
+ unit='deg')
+
+ return generate
+
+@pytest.fixture()
+def random_lonlat():
+ def generate(size):
+ return np.random.uniform(0, 360, size) * u.deg, \
+ np.random.uniform(-89, 89, size) * u.deg
+
+ return generate
+
+@pytest.fixture()
+def mocs():
+ moc1 = {'1': [0]}
+ moc1_increased = {'1': [0, 1, 2, 17, 22]}
+ moc2 = {'1': [30]}
+ moc2_increased = {'1': [8, 30, 31, 28, 43]}
+
+ return dict(moc1=MOCSkyRegion.from_json(moc1),
+ moc1_increased=MOCSkyRegion.from_json(moc1_increased),
+ moc2=MOCSkyRegion.from_json(moc2),
+ moc2_increased=MOCSkyRegion.from_json(moc2_increased))
+
+@pytest.fixture()
+def mocs_op():
+ moc1 = MOCSkyRegion.from_json({
+ '0': [0, 2, 3, 4, 5]
+ })
+ moc2 = MOCSkyRegion.from_json({
+ '0': [0, 1, 7, 4, 3]
+ })
+ moc3 = MOCSkyRegion.from_json({
+ '0': [6, 3]
+ })
+ return dict(first=moc1, second=moc2, third=moc3)
+
+@pytest.fixture()
+def wcs():
+ config = dict(crpix=(500, 500), crval=(0, 0), cdelt=(-0.1, 0.1), shape=(1000, 1000))
+ dataset = make_example_dataset(config=config)
+ return dataset.wcs
+
+class TestMOC(object):
+ """Test involving the manipulation of MOC objects"""
+
+ def setup_class(self):
+ # Load from fits file
+ filename_galex = get_pkg_data_filename('shapes/tests/data/P-GALEXGR6-AIS-FUV.fits', package='regions')
+ self.galex = MOCSkyRegion.from_fits(filename_galex)
+ filename_sdss = get_pkg_data_filename('shapes/tests/data/P-SDSS9-r.fits', package='regions')
+ self.sdss = MOCSkyRegion.from_fits(filename_sdss)
+
+ @pytest.mark.parametrize("size", [
+ 1000, 10000, 50000
+ ])
+ def test_from_skycoord(self, random_skycoords, size):
+ skycoord = random_skycoords(size)
+ moc = MOCSkyRegion.from_skycoord(skycoord, max_depth=7)
+
+ @pytest.mark.parametrize("size", [
+ 1000, 10000, 50000
+ ])
+ def test_from_lonlat(self, random_lonlat, size):
+ lon, lat = random_lonlat(size)
+ moc = MOCSkyRegion.from_lonlat(lon=lon, lat=lat, max_depth=7)
+
+ def test_from_cells(self):
+ num_ipix = 1000
+ cells = np.zeros(num_ipix, dtype={
+ 'names':('ipix', 'depth', 'fully_covered'),
+ 'formats':(np.uint64, np.uint32, np.uint8),
+ })
+ depth = 12
+ npix = 12 * (4 ** depth)
+
+ for i in range(num_ipix):
+ cells["ipix"][i] = np.random.randint(npix, size=1)
+ cells["depth"][i] = depth
+ cells["fully_covered"] = 1
+
+ MOCSkyRegion.from_cells(cells)
+
+ def test_from_fits(self):
+ assert self.galex
+
+ def test_moc_serialize_and_from_json(self):
+ ipix_d = self.galex.serialize(format="json")
+ moc = MOCSkyRegion.from_json(ipix_d)
+ assert self.galex == moc
+
+ def test_write_to_json(self):
+ data = self.galex.write("to_json.txt", format='json')
+ import json
+ with open('to_json.txt', 'r') as f_in:
+ data = json.load(f_in)
+ moc_from_file = MOCSkyRegion.from_json(data)
+ assert self.galex == moc_from_file
+
+ def test_write_to_fits(self):
+ data = self.galex.write("galex.fits", format='fits')
+ moc_from_file = MOCSkyRegion.from_fits("galex.fits")
+ assert self.galex == moc_from_file
+
+ def test_serialize_to_fits(self):
+ hdulist = self.galex.serialize(format='fits')
+ assert isinstance(hdulist, fits.hdu.hdulist.HDUList)
+
+ def test_serialize_to_json(self):
+ data = self.galex.serialize(format='json')
+ assert isinstance(data, dict)
+
+ def test_contains(self):
+ level = 4
+ size = 20
+ ipix = np.random.randint(0, 12*4**level, size)
+ npix = np.arange(12*4**level)
+ ipix_out = np.setdiff1d(npix, ipix)
+
+ moc = MOCSkyRegion.from_json({str(level): list(ipix)})
+
+ hp = HEALPix(nside=(1 << level), order='nested', frame=ICRS())
+ lon_in, lat_in = hp.healpix_to_lonlat(ipix)
+ lon_out, lat_out = hp.healpix_to_lonlat(ipix_out)
+
+ test_in = moc.contains(ra=lon_in, dec=lat_in)
+ assert test_in.all()
+ test_out = moc.contains(ra=lon_out, dec=lat_out)
+ assert not test_out.any()
+
+ def test_add_neighbours(self, mocs):
+ assert mocs['moc1'].add_neighbours() == mocs['moc1_increased']
+ assert mocs['moc2'].add_neighbours() == mocs['moc2_increased']
+
+ def test_remove_neighbours(self, mocs):
+ assert mocs['moc1_increased'].remove_neighbours() == mocs['moc1']
+ assert mocs['moc2_increased'].remove_neighbours() == mocs['moc2']
+
+ def test_neighbours(self, mocs):
+ moc1 = copy.deepcopy(mocs['moc1'])
+ moc2 = copy.deepcopy(mocs['moc2'])
+ assert moc1.add_neighbours().remove_neighbours() == mocs['moc1']
+ assert moc2.add_neighbours().remove_neighbours() == mocs['moc2']
+
+ def test_sky_fraction(self):
+ moc = MOCSkyRegion.from_json({
+ '0': [0, 1, 2, 3, 4, 5]
+ })
+ assert moc.sky_fraction == 0.5
+
+ def test_union(self, mocs_op):
+ assert mocs_op['first'].union(mocs_op['second'], mocs_op['third']) == MOCSkyRegion.from_json({
+ '0': [0, 1, 2, 3, 4, 5, 6, 7]
+ })
+
+ def test_intersection(self, mocs_op):
+ assert mocs_op['first'].intersection(mocs_op['second'], mocs_op['third']) == MOCSkyRegion.from_json({
+ '0': [3]
+ })
+
+ def test_difference(self, mocs_op):
+ assert mocs_op['first'].difference(mocs_op['second'], mocs_op['third']) == MOCSkyRegion.from_json({
+ '0': [2, 5]
+ })
+
+ def test_degrade_moc(self):
+ precise_moc = MOCSkyRegion.from_json({
+ '1': [4, 21]
+ })
+ degraded_moc = precise_moc.degrade_to_depth(0)
+ assert degraded_moc == MOCSkyRegion.from_json({'0': [1, 5]})
+
+ def test_complement(self):
+ assert self.galex.complement().complement() == self.galex
+
+ @pytest.mark.parametrize("pos_x, pos_y, expected_result", [
+ (500, 500, False),
+ (700, 500, False),
+ (0, 0, True),
+ (1000, 1000, True),
+ ])
+ def test_contains_pix(self, wcs, pos_x, pos_y, expected_result):
+ pixel = PixCoord(pos_x, pos_y)
+
+ reg = self.galex.to_pixel(wcs)
+ assert reg.contains(pixel) == expected_result
+
+ def test_to_sky(self, wcs):
+ reg_px = self.galex.to_pixel(wcs)
+ reg_cast_back_sky = reg_px.to_sky(wcs)
+
+ assert reg_cast_back_sky == self.galex
+
+ @pytest.mark.skipif('not HAS_MATPLOTLIB')
+ def test_as_artist_moc(self, wcs):
+ import matplotlib.pyplot as plt
+ fig = plt.figure(111, figsize=(10, 10))
+
+ ax = fig.add_subplot(1, 1, 1, projection=wcs)
+ reg = self.galex.to_pixel(wcs)
+ reg.plot(ax=ax, alpha=0.5, fill=True, linewidth=0.5, color='r')
+
+ plt.xlabel('ra')
+ plt.ylabel('dec')
+ plt.title('Coverage of GALEX')
+ plt.grid(color="black", linestyle="dotted")
+ plt.clf()
+
+ @pytest.mark.skipif('not HAS_MATPLOTLIB')
+ def test_plot(self):
+ import matplotlib.pyplot as plt
+ fig = plt.figure(111, figsize=(10, 10))
+
+ with WCS(fig,
+ fov=50 * u.deg,
+ center=SkyCoord(0, 90, unit='deg', frame='icrs'),
+ coordsys="icrs",
+ rotation=Angle(0, u.degree),
+ projection="AIT") as wcs:
+ ax = fig.add_subplot(1, 1, 1, projection=wcs)
+ self.galex.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, linewidth=0, color='r')
+ self.galex.border(ax=ax, wcs=wcs, color='r')
+
+ plt.xlabel('ra')
+ plt.ylabel('dec')
+ plt.title('Coverage of GALEX')
+ plt.grid(color="black", linestyle="dotted")
+ plt.clf()
+
+ def test_boundaries(self):
+ moc = self.galex.degrade_to_depth(6)
+ boundaries_l = moc.get_boundaries()
+
+ def test_to_mask(self):
+ config = dict(crpix=(50, 50), crval=(0, 0), cdelt=(-0.1, 0.1), shape=(100, 100))
+ dataset = make_example_dataset(config=config)
+
+ reg_px = self.galex.to_pixel(dataset.wcs)
+ mask = reg_px.to_mask()
diff --git a/regions/shapes/tests/test_moc_from_polygon.py b/regions/shapes/tests/test_moc_from_polygon.py
new file mode 100644
index 00000000..fa464401
--- /dev/null
+++ b/regions/shapes/tests/test_moc_from_polygon.py
@@ -0,0 +1,161 @@
+from ..moc import MOCSkyRegion, WCS
+
+from astropy import units as u
+from astropy.coordinates import SkyCoord
+
+import numpy as np
+
+def test_create_from_polygon():
+ lon = [83.71315909, 83.71378887, 83.71297292, 83.71233919] * u.deg
+ lat = [-5.44217436,-5.44298864, -5.44361751, -5.4428033] * u.deg
+ inside = SkyCoord(-15, -15, unit="deg", frame="icrs")
+
+ moc = MOCSkyRegion.from_polygon(lon=lon, lat=lat, max_depth=20, inside=None)
+
+ truth_ipix_d = {'17': [89921647231],
+ '18': [359686588915,
+ 359686588918,
+ 359686588919,
+ 359686588921,
+ 359686588923,
+ 359686589268,
+ 359686589269,
+ 359686589608,
+ 359686589610,
+ 359686589952],
+ '19': [1438746355657,
+ 1438746355659,
+ 1438746357060,
+ 1438746357061,
+ 1438746357063,
+ 1438746358408,
+ 1438746358410,
+ 1438746359814],
+ '20': [5754985422603,
+ 5754985422606,
+ 5754985422607,
+ 5754985422618,
+ 5754985422619,
+ 5754985422622,
+ 5754985422623,
+ 5754985422625,
+ 5754985422627,
+ 5754985422633,
+ 5754985422635,
+ 5754985422664,
+ 5754985422665,
+ 5754985422666,
+ 5754985422667,
+ 5754985422668,
+ 5754985422669,
+ 5754985422670,
+ 5754985422671,
+ 5754985422680,
+ 5754985422681,
+ 5754985422682,
+ 5754985422683,
+ 5754985422684,
+ 5754985422685,
+ 5754985422686,
+ 5754985422687,
+ 5754985422721,
+ 5754985422724,
+ 5754985422725,
+ 5754985422726,
+ 5754985422727,
+ 5754985422732,
+ 5754985422733,
+ 5754985422734,
+ 5754985422735,
+ 5754985422756,
+ 5754985422757,
+ 5754985422758,
+ 5754985422759,
+ 5754985422765,
+ 5754985422767,
+ 5754985428229,
+ 5754985428231,
+ 5754985428237,
+ 5754985428248,
+ 5754985428249,
+ 5754985428250,
+ 5754985428251,
+ 5754985428272,
+ 5754985428273,
+ 5754985428274,
+ 5754985428275,
+ 5754985428276,
+ 5754985428277,
+ 5754985428278,
+ 5754985428279,
+ 5754985428320,
+ 5754985428321,
+ 5754985428322,
+ 5754985428323,
+ 5754985428324,
+ 5754985428325,
+ 5754985428326,
+ 5754985428327,
+ 5754985428336,
+ 5754985428337,
+ 5754985428338,
+ 5754985428339,
+ 5754985428340,
+ 5754985428341,
+ 5754985428342,
+ 5754985428343,
+ 5754985433608,
+ 5754985433609,
+ 5754985433610,
+ 5754985433611,
+ 5754985433612,
+ 5754985433613,
+ 5754985433614,
+ 5754985433615,
+ 5754985433636,
+ 5754985433637,
+ 5754985433638,
+ 5754985433639,
+ 5754985433644,
+ 5754985433645,
+ 5754985433646,
+ 5754985433647,
+ 5754985433656,
+ 5754985433658,
+ 5754985433744,
+ 5754985433746,
+ 5754985433752,
+ 5754985433754,
+ 5754985433755,
+ 5754985433776,
+ 5754985433777,
+ 5754985433778,
+ 5754985433779,
+ 5754985433784,
+ 5754985433785,
+ 5754985433786,
+ 5754985433787,
+ 5754985439248,
+ 5754985439249,
+ 5754985439250,
+ 5754985439251,
+ 5754985439252,
+ 5754985439254,
+ 5754985439260,
+ 5754985439262,
+ 5754985439264,
+ 5754985439265,
+ 5754985439266,
+ 5754985439267,
+ 5754985439268,
+ 5754985439269,
+ 5754985439270,
+ 5754985439271,
+ 5754985439280,
+ 5754985439281,
+ 5754985439282,
+ 5754985439283,
+ 5754985439284]}
+
+ moc_truth = MOCSkyRegion.from_json(truth_ipix_d)
+ assert(moc == moc_truth)
diff --git a/setup.py b/setup.py
index 4326ec44..d2e4c749 100755
--- a/setup.py
+++ b/setup.py
@@ -107,7 +107,7 @@
description=DESCRIPTION,
scripts=scripts,
setup_requires=['numpy'],
- install_requires=['numpy', 'astropy', 'astropy-healpix', 'six'],
+ install_requires=['numpy', 'astropy', 'astropy-healpix', 'six', 'networkx', 'spherical-geometry'],
extras_require=dict(
plot=[
'matplotlib',