Swith to cdshealpix for HEALPix polygon query

[lpsinger]

It seems to be more accurate than Healpy.

Here's an illustration. The example script plot_filling_factor.py calculates the filling factor of a grid of rectangular survey tiles as the density of the tiling increases. One expects the filling factor to vary monotonically with the density of the grid. However, this is not the case using hp.query_polygon:

However, with this patch to switch to cdshealpix.nested.polygon_search, it is monotonic, or very nearly:

[lpsinger]

@zonca, do you know why the cdshealpix polygon query might be more or less accurate than the healpy one?

[codecov[bot]]

Codecov Report

Merging #49 (594dd1f) into main (1e29f2f) will increase coverage by 0.35%. The diff coverage is 83.33%.

@@            Coverage Diff             @@
##             main      #49      +/-   ##
==========================================
+ Coverage   26.74%   27.09%   +0.35%     
==========================================
  Files          28       28              
  Lines        1350     1351       +1     
==========================================
+ Hits          361      366       +5     
+ Misses        989      985       -4     

Impacted Files	Coverage Δ
dorado/scheduling/fov.py	94.44% <83.33%> (+6.56%)	:arrow_up:
dorado/scheduling/scripts/simsurvey.py	9.09% <0.00%> (+0.08%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 1e29f2f...594dd1f. Read the comment docs.

[zonca]

thanks @lpsinger, @mreineck would you like to look into this issue?

[mreineck]

I'm not sure about all the details but you may be comparing apples and oranges here...

If I read the CDS docs correctly, their polygon search returns all pixels that overlap with the polygon. The polygon query of healpy by default only returns pixels whose center lies within the polygon, This can be changed by setting the inclusive parameter to True.

Please excuse me if I'm on the wrong track ...

[lpsinger]

@mreineck, I tried setting inclusive=True, and it did not make a difference.

[mreineck]

If you don't see any difference with inclusive=True, I expect that there is some bug in the calling chain to Healpix C++ where this setting is lost. I'll try to do some digging, but it may take time.

[mreineck]

It would be great if we could look at a map highlighting the pixels that are in the CDS query_polygon results, but not in the healpy ones. Unfortunately I'm having trouble installing all the required dependencies locally.

[lpsinger]

Unfortunately I'm having trouble installing all the required dependencies locally.

Which dependencies are causing you trouble?

[mreineck]

CPLEX originally, but I just commented it in pyproject.toml; it's apparently not needed for the script.

[lpsinger]

CPLEX originally, but I just commented it in pyproject.toml; it's apparently not needed for the script.

Interesting, what was the error message? CPLEX doesn't support Python 3.9 yet, FYI. You'd need to use 3.7 or 3.8 for this project.

[mreineck]

I'm on vanilla Debian testing which comes with Python 3.9 and I don't use virtual environments (yes, I know I'm strange :P). Also, a superficially glance at the CPLEX site seemed to suggest that I need an account to get the code, so my first choice was to avoid this dependency.

[lpsinger]

I'm on vanilla Debian testing which comes with Python 3.9 and I don't use virtual environments (yes, I know I'm strange :P).

Not that strange. That's a very solid Python runtime, in my opinion.

Also, a superficially glance at the CPLEX site seemed to suggest that I need an account to get the code, so my first choice was to avoid this dependency.

Yes, we have one commercial dependency. You don't need an account to install it.

Anyway, here's a more self-contained example with fewer dependencies. You can see that the Healpy version with inclusive=True has a few false positive pixels compared to the cdshealpix version.

# Adapted from https://github.com/nasa/dorado-scheduling/blob/main/dorado/scheduling/fov.py
#
# Copyright © 2020 United States Government as represented by the Administrator
# of the National Aeronautics and Space Administration. No copyright is claimed
# in the United States under Title 17, U.S. Code. All Other Rights Reserved.
#
# SPDX-License-Identifier: NASA-1.3
#
from astropy.coordinates import ICRS, SkyCoord, UnitSphericalRepresentation
from astropy_healpix import HEALPix
from astropy import units as u
import healpy as hp
from cdshealpix.nested import polygon_search
import ligo.skymap.plot
from matplotlib import pyplot as plt
import numpy as np

__all__ = ('FOV',)

class FOV:

    def __init__(self, representation):
        self._representation = representation

    @classmethod
    def from_rectangle(cls, width, height=None):
        """Create a rectangular field of view.

        Parameters
        ----------
        width : :class:`astropy.units.Quantity`
            The width of the rectangle.
        height : :class:`astropy.units.Quantity`
            The height of the rectangle. If omitted, then the field of view is
            a square.

        Returns
        -------
        :class:`FOV`

        """
        if height is None:
            height = width
        lon = np.asarray([0.5, 0.5, -0.5, -0.5]) * width
        lat = np.asarray([0.5, -0.5, -0.5, 0.5]) * height
        return cls(UnitSphericalRepresentation(lon, lat))

    def footprint(self, center=SkyCoord(0*u.deg, 0*u.deg), roll=0*u.deg):
        """Get the footprint of the FOV at a given orientation.

        Parameters
        ----------
        center : :class:`astropy.coordinates.SkyCoord`
            The center of the field of view. If omitted, the default is
            RA=0°, Dec=0°.
        roll : :class:`astropy.coordinates.SkyCoord`
            The roll of the field of view. If omitted, the default is 0°.

        Returns
        -------
        vertices : `astropy.coordinates.SkyCoord`
            The coordinates of the vertices of the field of view.

        Examples
        --------
        First, some imports:

        >>> from astropy.coordinates import ICRS, SkyCoord
        >>> from astropy import units as u
        >>> from astropy_healpix import HEALPix

        Now, get the footprint for the default pointing:

        >>> fov = FOV.from_rectangle(50 * u.deg)
        >>> fov.footprint().icrs
        <SkyCoord (ICRS): (ra, dec) in deg
            [( 25.,  25.), ( 25., -25.), (335., -25.), (335.,  25.)]>

        Get the footprint for a spcific pointing:

        >>> fov.footprint(SkyCoord(0*u.deg, 20*u.deg), 15*u.deg).icrs
        <SkyCoord (ICRS): (ra, dec) in deg
            [( 35.73851095,  34.84634609), ( 15.41057822, -11.29269295),
             (331.35544934,  -0.54495686), (336.46564791,  49.26075815)]>

        """
        frame = center.skyoffset_frame(roll)[..., np.newaxis]
        representation = self._representation

        # FIXME: Astropy does not automatically broadcast frame attributes.
        # Remove once https://github.com/astropy/astropy/issues/8812 is fixed.
        shape = np.broadcast_shapes(representation.shape, frame.shape)
        frame = np.broadcast_to(frame, shape, subok=True)
        representation = np.broadcast_to(self._representation, shape,
                                         subok=True)

        # FIXME: Astropy does not support realizing frames with representations
        # that are broadcast views. Remove once
        # https://github.com/astropy/astropy/issues/11572 is fixed.
        representation = representation.copy()

        return SkyCoord(frame.realize_frame(representation))

    @staticmethod
    def _polygon_search_internal(healpix, vertices):
        ipix, _, _ = polygon_search(
            vertices.lon, vertices.lat, depth=healpix.level, flat=True)
        ipix = ipix.astype(np.intp)
        if healpix.order == 'ring':
            ipix = healpix.nested_to_ring(ipix)
        return ipix

    def footprint_healpix_cds(self, healpix, *args, **kwargs):
        """Get the HEALPix footprint at a given orientation.

        Parameters
        ----------
        center : :class:`astropy.coordinates.SkyCoord`
            The center of the field of view.
        rotate : class:`astropy.units.Quantity`
            The position angle (optional, default 0 degrees).

        Returns
        -------
        :class:`np.ndarray`
            An array of HEALPix indices contained within the footprint.

        """
        vertices = self.footprint(*args, **kwargs).transform_to(
            healpix.frame).represent_as(UnitSphericalRepresentation)
        return self._polygon_search_internal(healpix, vertices)

    def footprint_healpix_healpy(self, healpix, *args, **kwargs):
        vertices = self.footprint(*args, **kwargs)
        xyz = vertices.transform_to(healpix.frame).cartesian.xyz.value
        return hp.query_polygon(healpix.nside, xyz.T,
                                nest=(healpix.order == 'nested'))

    def footprint_healpix_healpy_inclusive(self, healpix, *args, **kwargs):
        vertices = self.footprint(*args, **kwargs)
        xyz = vertices.transform_to(healpix.frame).cartesian.xyz.value
        return hp.query_polygon(healpix.nside, xyz.T, inclusive=True,
                                nest=(healpix.order == 'nested'))

fov = FOV.from_rectangle(6.8 * u.deg)
healpix = HEALPix(nside=32, frame=ICRS())

methods = [fov.footprint_healpix_cds,
           fov.footprint_healpix_healpy,
           fov.footprint_healpix_healpy_inclusive]

fig, axs = plt.subplots(
    3, 1, figsize=(10, 30),
    subplot_kw=dict(projection='astro zoom', radius=6 * u.deg))
for ax, method in zip (axs, methods):
    ax.set_title(method.__name__)
    coords = fov.footprint().icrs
    ax.add_patch(plt.Polygon(
        np.column_stack((coords.ra.deg, coords.dec.deg)),
        transform=ax.get_transform('world'),
        facecolor='none', edgecolor='blue'))
    ipix = method(healpix)
    for boundaries in healpix.boundaries_skycoord(ipix, 1):
        ax.add_patch(plt.Polygon(
            np.column_stack((boundaries.ra.deg, boundaries.dec.deg)),
            transform=ax.get_transform('world'),
            facecolor='none', edgecolor='red'))

[mreineck]

Thanks a lot, that script works out of the box!

It is possible to avoid the false positives by passing a higher fact parameter to query_polygon (the default is 4, and it must be a power of 2). In this concrete example, I get the same result as cdshealpix already for fact=8.

However the question remains why the polygon query behaves so differently in your filling factor tests ... since the inclusive search in healpy produces false positives, I would have expected that you actually get even higher fill factors than cdshealpix when you switch inclusive searches on. This is really strange.

[lpsinger]

It is possible to avoid the false positives by passing a higher fact parameter to query_polygon (the default is 4, and it must be a power of 2). In this concrete example, I get the same result as cdshealpix already for fact=8.

However the question remains why the polygon query behaves so differently in your filling factor tests ... since the inclusive search in healpy produces false positives, I would have expected that you actually get even higher fill factors than cdshealpix when you switch inclusive searches on. This is really strange.

This is what the result looks like using Healpy and inclusive=True, fact=8. Not much better.