Feat add box and zone

.github/workflows/test.yml

@@ -20,6 +20,7 @@ jobs:
     # See https://github.com/pypa/manylinux for this particular container:
     # * CPython 3.8, 3.9, ... installed in /opt/python/<python tag>-<abi tag>
     container: quay.io/pypa/manylinux2014_x86_64
+    env: {ACTIONS_ALLOW_USE_UNSECURE_NODE_VERSION: true} # Allow using Node16 actions required for CentOS7
     # We are now in CentOS 7 64 bits
     steps:
      - name: "Checkout the full project"

.gitignore

@@ -9,3 +9,4 @@ docs/_build
 docs/_templates
 publish_docs.sh
 .idea
+*.so

.pre-commit-config.yaml

@@ -35,3 +35,11 @@ repos:
     hooks:
      - id: ruff
        fail_fast: true
+  - repo: local
+    hooks:
+      - id: cargofmt
+        name: cargofmt
+        entry: rustfmt
+        files: \.rs$
+        language: system
+

CHANGELOG.md

@@ -2,9 +2,21 @@
 
 ## Unreleased
 
-### Fix
+### Added
+
+* new feature `box_search`
+* `cdshaelpix.nested.vertices` can now take depth as a `numpy.ndarray` instead of only
+  accepting a single depth
+
+### Fixed
 
-* bilinear_interpolation now accepts longitudes and latitudes with nan values (will be a masked value in the output) [#22]
+* :warning: breaking change. All `fully_covered` returned values have changed from arrays
+  of 0 and 255 to boolean arrays. This concerns: `box_search`, `cone_search`,
+  `polygon_search`, `elliptical_cone_search`, and `elliptic-search`.
+* bilinear_interpolation now accepts longitudes and latitudes with nan values
+  (will be a masked value in the output) [#22]
+* `nested.healpix_to_lonlat`, failed into rust panic for `dx=1` or `dy=1`. This is
+  now indicated in the documentation and is catched in a `ValueError` on the python side.
 
 ## 0.6.5
 

docs/conf.py

@@ -230,10 +230,10 @@
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    "python": ("https://docs.python.org/", None),
-    "astropy": ("http://docs.astropy.org/en/latest/", None),
-    "numpy": ("http://docs.scipy.org/doc/numpy/", None),
-    "matplotlib": ("https://matplotlib.org/", None),
+    "python": ("https://docs.python.org/3/", None),
+    "astropy": ("https://docs.astropy.org/en/latest/", None),
+    "numpy": ("https://numpy.org/doc/stable/", None),
+    "matplotlib": ("https://matplotlib.org/stable/", None),
 }
 
 

docs/examples/examples.rst

@@ -1,26 +1,54 @@
+********
 Examples
-========
+********
+
+Getting coordinates in a HEALPix cell
+=====================================
+
+In this example, we chose a HEALPix cell, and we plot its center
+
+.. plot:: examples/healpix_to_lonlat.py
+    :include-source:
+
+
+Getting the healpix cells contained in a region of the sky
+==========================================================
 
 Cone search
 -----------
 
-.. plot:: examples/cone_search.py
+.. plot:: examples/search_methods/cone_search.py
     :include-source:
 
 Elliptical cone search
 ----------------------
 
-.. plot:: examples/elliptic_search.py
+.. plot:: examples/search_methods/elliptic_search.py
     :include-source:
 
 Polygon search
 --------------
 
-.. plot:: examples/polygon_search.py
+.. plot:: examples/search_methods/polygon_search.py
+    :include-source:
+
+Box search
+----------
+
+.. plot:: examples/search_methods/box_search.py
+    :include-source:
+
+Zone search
+-----------
+
+In this example, we get the ``ipix`` and ``depth`` in a zone and plot them by combining
+:external:ref:`~cdshealpix.nested.vertices`_ with :external:ref:`~matplotlib.path.Polygon`_
+
+.. plot:: examples/search_methods/zone_search.py
     :include-source:
 
 Notebook examples
------------------
+=================
 
 .. nbgallery::
     ../_collections/notebooks/coordinate_conversion.ipynb

docs/examples/healpix_to_lonlat.py

@@ -0,0 +1,28 @@
+"""Plot a cell and its center."""
+
+from cdshealpix.nested import healpix_to_lonlat, vertices
+
+import astropy.units as u
+import matplotlib.pyplot as plt
+from matplotlib.patches import Polygon
+import numpy as np
+
+# choose a cell
+ipix = 5
+depth = 0
+
+fig = plt.figure()
+ax = fig.add_subplot(projection="aitoff")
+ax.grid(visible=True)
+# plot the center of the cell
+lon, lat = healpix_to_lonlat(ipix, depth)
+ax.scatter(lon, lat, color="purple")
+# plot the cell
+border_lon, border_lat = vertices(ipix, depth, step=5)
+border_lon = border_lon.wrap_at(180 * u.deg).radian[0]
+border_lat = border_lat.radian[0]
+polygon = Polygon(
+    np.column_stack((border_lon, border_lat)), fill=False, edgecolor="green", hatch="\\"
+)
+ax.add_patch(polygon)
+plt.show()

docs/examples/search_methods/box_search.py

@@ -0,0 +1,60 @@
+"""Demonstration of the cone_search function."""
+
+# Astropy tools
+import astropy.units as u
+
+# For plots
+import matplotlib.pyplot as plt
+from astropy.coordinates import Angle, Latitude, Longitude, SkyCoord
+
+# Moc and HEALPix tools
+from cdshealpix.nested import box_search
+from mocpy import MOC, WCS as mocpy_WCS  # noqa: N811
+
+max_depth = 10
+
+ipix, depth, fully_covered = box_search(
+    lon=Longitude(0, u.deg),
+    lat=Latitude(0, u.deg),
+    a=2 * u.deg,
+    b=1 * u.deg,
+    angle=0 * u.deg,
+    depth=max_depth,
+)
+
+
+moc_fully_covered = MOC.from_healpix_cells(
+    ipix[fully_covered], depth[fully_covered], max_depth
+)
+moc_partially_covered = MOC.from_healpix_cells(
+    ipix[~fully_covered], depth[~fully_covered], max_depth
+)
+
+# plot
+fig = plt.figure(111, figsize=(10, 10))
+# Define a astropy WCS from the mocpy.WCS class
+with mocpy_WCS(
+    fig,
+    fov=5 * u.deg,
+    center=SkyCoord(0, 0, unit="deg", frame="icrs"),
+    coordsys="icrs",
+    rotation=Angle(0, u.degree),
+    projection="AIT",
+) as wcs:
+    ax = fig.add_subplot(1, 1, 1, projection=wcs)
+    # Call fill with a matplotlib axe and the `~astropy.wcs.WCS` wcs object.
+    moc_fully_covered.fill(
+        ax=ax, wcs=wcs, alpha=0.5, fill=True, color="green", label="fully inside"
+    )
+    moc_partially_covered.fill(
+        ax=ax, wcs=wcs, alpha=0.5, fill=True, color="red", label="crosses the border"
+    )
+
+plt.xlabel("ra")
+plt.ylabel("dec")
+plt.legend()
+plt.title(
+    "Box search, accessing the cells that are partially covered by the sky region"
+)
+plt.grid(color="black", linestyle="dotted")
+plt.show()

docs/examples/cone_search.py → docs/examples/search_methods/cone_search.py

@@ -8,7 +8,7 @@
 from astropy.coordinates import Angle, Latitude, Longitude, SkyCoord
 
 # Moc and HEALPix tools
-from cdshealpix import cone_search
+from cdshealpix.nested import cone_search
 from mocpy import MOC, WCS
 
 max_depth = 10

docs/examples/search_methods/zone_search.py

@@ -0,0 +1,35 @@
+"""Get the HEALPix cells in a zone, and plot all of them as matplotlib polygons."""
+
+from cdshealpix.nested import zone_search, vertices
+
+from astropy.coordinates import Longitude, Latitude
+
+import astropy.units as u
+import numpy as np
+from matplotlib.collections import PolyCollection
+import matplotlib.pyplot as plt
+
+lon_min = Longitude(-10 * u.deg)
+lat_min = Latitude(10 * u.deg)
+lon_max = Longitude(60 * u.deg)
+lat_max = Latitude(60 * u.deg)
+
+ipix, depth, _ = zone_search(lon_min, lat_min, lon_max, lat_max, depth=4)
+
+n_cells = len(ipix)
+
+lons, lats = vertices(ipix, depth, step=4)
+
+lons = lons.wrap_at(180 * u.deg).radian
+lats = lats.radian
+
+path_polygons = np.array([np.column_stack((lon, lat)) for lon, lat in zip(lons, lats)])
+
+polygons = PolyCollection(path_polygons)
+
+fig = plt.figure()
+ax = fig.add_subplot(projection="aitoff")
+ax.grid(visible=True)
+ax.add_collection(polygons)
+
+plt.show()

docs/index.rst

@@ -7,7 +7,7 @@ Welcome to cdshealpix's documentation!
 ======================================
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 3
    :caption: Contents:
 
    install

pyproject.toml

@@ -65,3 +65,8 @@ max-string-length = 20
 
 [tool.ruff.pydocstyle]
 convention = "numpy"  # Accepts: "google", "numpy", or "pep257"
+
+[tool.pytest.ini_options]
+addopts = "--doctest-modules"
+doctest_optionflags = "NORMALIZE_WHITESPACE"
+testpaths = "python"