wwt_data_formats/place.py: add support for constellations

Here we import WWT's constellation database and emulate its logic for
figuring out the constellation associated with a given RA/Dec. This will
certainly disagree with Astropy's implementation for corner cases, but
for WWT purposes we should emulate WWT logic.

MANIFEST.in

@@ -6,6 +6,7 @@ include *.yml
 graft docs
 prune docs/_build
 
+graft wwt_data_formats/data
 graft wwt_data_formats/tests
 
 global-exclude *~

docs/api/wwt_data_formats.place.Place.rst

@@ -61,6 +61,7 @@ Place
       ~Place.notify_change
       ~Place.observe
       ~Place.on_trait_change
+      ~Place.set_ra_dec
       ~Place.set_trait
       ~Place.setup_instance
       ~Place.to_xml
@@ -74,6 +75,7 @@ Place
       ~Place.traits
       ~Place.unobserve
       ~Place.unobserve_all
+      ~Place.update_constellation
       ~Place.write_xml
 
    .. rubric:: Attributes Documentation
@@ -127,6 +129,7 @@ Place
    .. automethod:: notify_change
    .. automethod:: observe
    .. automethod:: on_trait_change
+   .. automethod:: set_ra_dec
    .. automethod:: set_trait
    .. automethod:: setup_instance
    .. automethod:: to_xml
@@ -140,4 +143,5 @@ Place
    .. automethod:: traits
    .. automethod:: unobserve
    .. automethod:: unobserve_all
+   .. automethod:: update_constellation
    .. automethod:: write_xml

wwt_data_formats/imageset.py

@@ -548,12 +548,10 @@ class ImageSet(LockedXmlTraits, UrlContainer):
                 self.rotation_deg = -self.rotation_deg
 
         if place is not None:
-            place.data_set_type = DataSetType.SKY
+            place.set_ra_dec(center_ra_deg / 15.0, center_dec_deg)
             place.rotation_deg = (
                 0.0  # I think this is better than propagating the image rotation?
             )
-            place.ra_hr = center_ra_deg / 15.0
-            place.dec_deg = center_dec_deg
             # It is hardcoded that in sky mode, zoom_level = height of client FOV * 6.
             place.zoom_level = height * scale_y * fov_factor * 6
 

wwt_data_formats/place.py

@@ -1,10 +1,11 @@
 # -*- mode: python; coding: utf-8 -*-
-# Copyright 2019-2020 the .NET Foundation
+# Copyright 2019-2022 the .NET Foundation
 # Licensed under the MIT License.
 
-"""A place that a WWT user can visit.
-
 """
+A place that a WWT user can visit.
+"""
+
 from __future__ import absolute_import, division, print_function
 
 __all__ = """
@@ -12,6 +13,7 @@ Place
 """.split()
 
 from argparse import Namespace
+from collections import namedtuple
 from traitlets import Float, Instance, Int, Unicode, UseEnum
 
 from . import LockedXmlTraits, XmlSer
@@ -31,9 +33,15 @@ class Place(LockedXmlTraits, UrlContainer):
     dec_deg = Float(0.0).tag(xml=XmlSer.attr("Dec"), xml_if_sky_type_is=True)
     latitude = Float(0.0).tag(xml=XmlSer.attr("Lat"), xml_if_sky_type_is=False)
     longitude = Float(0.0).tag(xml=XmlSer.attr("Lng"), xml_if_sky_type_is=False)
+
     constellation = UseEnum(Constellation, default_value=Constellation.UNSPECIFIED).tag(
         xml=XmlSer.attr("Constellation")
     )
+    """
+    The constellation associated with this place's sky position, if it has one.
+    Use :meth:`set_ra_dec` to compute this correctly (according to WWT).
+    """
+
     classification = UseEnum(
         Classification, default_value=Classification.UNSPECIFIED
     ).tag(xml=XmlSer.attr("Classification"))
@@ -130,3 +138,183 @@ class Place(LockedXmlTraits, UrlContainer):
         if self.foreground_image_set is not None:
             return self.foreground_image_set
         return self.image_set
+
+    def set_ra_dec(self, ra_hr, dec_deg):
+        """
+        Update the sky coordinates associated with this Place and associated
+        metadata.
+
+        Parameters
+        ----------
+        ra_hr : number
+            The right ascension, in hours
+        dec_deg : number
+            The declination, in degrees
+
+        Returns
+        -------
+        Self, for chaining.
+
+        Notes
+        -----
+        Beyond simply setting the associated fields, this method will set the
+        :attr:`data_set_type` to Sky, compute the correct setting for
+        :attr:`constellation`, and clear :attr:`latitude` and :attr:`longitude`.
+        The constellation computation is done by emulating WWT's internal
+        computation, which may give different results than other methods in
+        (literal) corner cases. In cases where precision matters and WWT
+        compatibility does not, :func:`astropy.coordinates.get_constellation`
+        should be preferred.
+        """
+        db = _get_iau_constellations()
+        self.ra_hr = ra_hr
+        self.dec_deg = dec_deg
+        self.data_set_type = DataSetType.SKY
+        self.latitude = 0
+        self.longitude = 0
+        self.constellation = db.find_constellation_for_point(ra_hr, dec_deg)
+        return self
+
+    def update_constellation(self):
+        """
+        Update the constellation associated with this Place to agree with its
+        current RA and declination, if set.
+
+        Returns
+        -------
+        Self, for chaining.
+
+        Notes
+        -----
+        If this object has a :attr:`data_set_type` of Sky, this method will
+        update :attr:`constellation` to be correct given its :attr:`ra_hr` and
+        :attr:`dec_deg`, and clear :attr:`latitude` and :attr:`longitude`.
+        Otherwise, the constellation, RA, and Dec will be cleared.
+
+        The constellation computation is done by emulating WWT's internal
+        computation, which may give different results than other methods in
+        (literal) corner cases. In cases where precision matters and WWT
+        compatibility does not, :func:`astropy.coordinates.get_constellation`
+        should be preferred.
+        """
+
+        if self.data_set_type == DataSetType.SKY:
+            db = _get_iau_constellations()
+            self.latitude = 0
+            self.longitude = 0
+            self.constellation = db.find_constellation_for_point(
+                self.ra_hr, self.dec_deg
+            )
+        else:
+            self.ra_hr = 0
+            self.dec_deg = 0
+            self.constellation = Constellation.UNSPECIFIED
+
+        return self
+
+
+Radec = namedtuple("Radec", "ra_hr dec_deg")
+
+
+class ConstellationDatabase(object):
+    _shapes = None
+
+    def __init__(self, data_stream):
+        # This implementation emulates the WWT `Constellations()` constructor
+        # with `boundry=true` and `noInterpollation=true`, which is what is used
+        # for the constellation-containment tests (`constellationCheck` AKA
+        # `Containment`). We have also edited the `constellations.txt` file,
+        # pulled from
+        # `http://www.worldwidetelescope.org/data/constellations.txt`, to remove
+        # the `I` (interpolated) entries, which increase the size of the data
+        # file substantially and are ignored in the check mode, and to convert
+        # leading spaces to trailing spaces, especially with separated sign
+        # characters.
+
+        prev_abbrev = ""
+        cur_points = None
+        shapes = {}
+
+        for line in data_stream:
+            ra_hr = float(line[:10])
+            dec_deg = float(line[11:22])
+            abbrev = line[23:27].strip()  # "SER1" and "SER2" are 4 letters
+
+            if abbrev != prev_abbrev:
+                cur_points = []
+                shapes[Constellation(abbrev)] = cur_points
+                prev_abbrev = abbrev
+                prev_ra = 0.0
+
+            if ra_hr - prev_ra > 12:
+                ra_hr -= 24
+            elif ra_hr - prev_ra < -12:
+                ra_hr += 24
+
+            cur_points.append(Radec(ra_hr, dec_deg))
+            prev_ra = ra_hr
+
+        self._shapes = shapes
+
+    def find_constellation_for_point(self, ra_hr, dec_deg):
+        """
+        Find the constellation corresponding to a given RA and Dec.
+
+        ra_hr : RA in hours
+        dec_deg : dec in degrees
+        """
+        # This implementation from WWT's `FindConstellationForPoint()` function.
+
+        if dec_deg > 88.402:
+            return Constellation.URSA_MINOR
+
+        for constellation, points in self._shapes.items():
+            inside = False
+            j = len(points) - 1
+
+            for i in range(len(points)):
+                c = points[i].dec_deg <= dec_deg and dec_deg < points[j].dec_deg
+                c = c or points[j].dec_deg <= dec_deg and dec_deg < points[i].dec_deg
+
+                if c:
+                    x = (points[j].ra_hr - points[i].ra_hr) * (
+                        dec_deg - points[i].dec_deg
+                    )
+                    x /= points[j].dec_deg - points[i].dec_deg
+                    if ra_hr < x + points[i].ra_hr:
+                        inside = not inside
+
+                j = i
+
+            if inside:
+                return constellation
+
+        if ra_hr > 0:
+            return self.find_constellation_for_point(ra_hr - 24, dec_deg)
+
+        # "Ursa Minor is tricky since it wraps around the poles. I[t] can evade
+        # the point in rect test."
+
+        if dec_deg > 65.5:
+            return Constellation.URSA_MINOR
+
+        if dec_deg < -65.5:
+            return Constellation.OCTANS
+
+        return Constellation.UNSPECIFIED
+
+
+_iau_constellation_data = None
+
+
+def _get_iau_constellations():
+    import os.path
+
+    global _iau_constellation_data
+
+    if _iau_constellation_data is None:
+        data = os.path.join(os.path.dirname(__file__), "data", "iau_constellations.txt")
+        with open(data, "rt") as f:
+            _iau_constellation_data = ConstellationDatabase(f)
+
+    return _iau_constellation_data

wwt_data_formats/tests/test_place.py

@@ -8,6 +8,7 @@ from xml.etree import ElementTree as etree
 
 from . import assert_xml_trees_equal
 from .. import imageset, place
+from ..enums import Constellation
 
 
 def test_basic_xml():
@@ -152,3 +153,20 @@ def test_nesting():
     pl.background_image_set.url = "http://example.com/background"
     observed_xml = pl.to_xml()
     assert_xml_trees_equal(expected_xml, observed_xml)
+
+
+def test_constellations():
+    SAMPLES = [
+        (23.99, 90, Constellation.URSA_MINOR),
+        (1.5, 82.5, Constellation.CEPHEUS),
+        (20.5, 41.5, Constellation.CYGNUS),
+        (17.0, -42.6, Constellation.SCORPIUS),
+        (0, -90, Constellation.OCTANS),
+        (6, -84.5, Constellation.MENSA),
+    ]
+
+    pl = place.Place()
+
+    for ra_hr, dec_deg, expected in SAMPLES:
+        pl.set_ra_dec(ra_hr, dec_deg)
+        assert pl.constellation == expected