Refactor of plotting contours

dependencies-examples.log

@@ -8,29 +8,29 @@ astropy==5.3.4
     # via
     #   healpy
     #   icecube-skyreader (setup.py)
-cachetools==5.3.1
+cachetools==5.3.2
     # via wipac-rest-tools
 certifi==2023.7.22
     # via requests
 cffi==1.16.0
     # via cryptography
-charset-normalizer==3.3.0
+charset-normalizer==3.3.2
     # via requests
-contourpy==1.1.1
+contourpy==1.2.0
     # via matplotlib
-cryptography==41.0.4
+cryptography==41.0.5
     # via pyjwt
 cycler==0.12.1
     # via matplotlib
-fonttools==4.43.1
+fonttools==4.44.0
     # via matplotlib
 healpy==1.16.6
     # via icecube-skyreader (setup.py)
 idna==3.4
     # via requests
 kiwisolver==1.4.5
     # via matplotlib
-matplotlib==3.8.0
+matplotlib==3.8.1
     # via
     #   healpy
     #   icecube-skyreader (setup.py)
@@ -50,7 +50,7 @@ packaging==23.2
     # via
     #   astropy
     #   matplotlib
-pandas==2.1.1
+pandas==2.1.2
     # via icecube-skyreader (setup.py)
 pillow==10.1.0
     # via matplotlib
@@ -99,7 +99,7 @@ urllib3==2.0.7
     # via
     #   requests
     #   wipac-rest-tools
-wipac-dev-tools==1.7.0
+wipac-dev-tools==1.7.1
     # via
     #   icecube-skyreader (setup.py)
     #   wipac-rest-tools

dependencies-mypy.log

@@ -8,21 +8,21 @@ astropy==5.3.4
     # via
     #   healpy
     #   icecube-skyreader (setup.py)
-cachetools==5.3.1
+cachetools==5.3.2
     # via wipac-rest-tools
 certifi==2023.7.22
     # via requests
 cffi==1.16.0
     # via cryptography
-charset-normalizer==3.3.0
+charset-normalizer==3.3.2
     # via requests
-contourpy==1.1.1
+contourpy==1.2.0
     # via matplotlib
-cryptography==41.0.4
+cryptography==41.0.5
     # via pyjwt
 cycler==0.12.1
     # via matplotlib
-fonttools==4.43.1
+fonttools==4.44.0
     # via matplotlib
 healpy==1.16.6
     # via icecube-skyreader (setup.py)
@@ -32,7 +32,7 @@ iniconfig==2.0.0
     # via pytest
 kiwisolver==1.4.5
     # via matplotlib
-matplotlib==3.8.0
+matplotlib==3.8.1
     # via
     #   healpy
     #   icecube-skyreader (setup.py)
@@ -53,7 +53,7 @@ packaging==23.2
     #   astropy
     #   matplotlib
     #   pytest
-pandas==2.1.1
+pandas==2.1.2
     # via icecube-skyreader (setup.py)
 pillow==10.1.0
     # via matplotlib
@@ -71,7 +71,7 @@ pyparsing==3.1.1
     # via matplotlib
 pypng==0.20220715.0
     # via qrcode
-pytest==7.4.2
+pytest==7.4.3
     # via
     #   icecube-skyreader (setup.py)
     #   pytest-mock
@@ -110,7 +110,7 @@ urllib3==2.0.7
     # via
     #   requests
     #   wipac-rest-tools
-wipac-dev-tools==1.7.0
+wipac-dev-tools==1.7.1
     # via
     #   icecube-skyreader (setup.py)
     #   wipac-rest-tools

dependencies-tests.log

@@ -10,13 +10,13 @@ astropy==5.3.4
     #   icecube-skyreader (setup.py)
 certifi==2023.7.22
     # via requests
-charset-normalizer==3.3.0
+charset-normalizer==3.3.2
     # via requests
-contourpy==1.1.1
+contourpy==1.2.0
     # via matplotlib
 cycler==0.12.1
     # via matplotlib
-fonttools==4.43.1
+fonttools==4.44.0
     # via matplotlib
 healpy==1.16.6
     # via icecube-skyreader (setup.py)
@@ -26,7 +26,7 @@ iniconfig==2.0.0
     # via pytest
 kiwisolver==1.4.5
     # via matplotlib
-matplotlib==3.8.0
+matplotlib==3.8.1
     # via
     #   healpy
     #   icecube-skyreader (setup.py)
@@ -47,7 +47,7 @@ packaging==23.2
     #   astropy
     #   matplotlib
     #   pytest
-pandas==2.1.1
+pandas==2.1.2
     # via icecube-skyreader (setup.py)
 pillow==10.1.0
     # via matplotlib
@@ -57,7 +57,7 @@ pyerfa==2.0.1.1
     # via astropy
 pyparsing==3.1.1
     # via matplotlib
-pytest==7.4.2
+pytest==7.4.3
     # via
     #   icecube-skyreader (setup.py)
     #   pytest-mock
@@ -83,5 +83,5 @@ tzdata==2023.3
     # via pandas
 urllib3==2.0.7
     # via requests
-wipac-dev-tools==1.7.0
+wipac-dev-tools==1.7.1
     # via icecube-skyreader (setup.py)

dependencies.log

@@ -10,21 +10,21 @@ astropy==5.3.4
     #   icecube-skyreader (setup.py)
 certifi==2023.7.22
     # via requests
-charset-normalizer==3.3.0
+charset-normalizer==3.3.2
     # via requests
-contourpy==1.1.1
+contourpy==1.2.0
     # via matplotlib
 cycler==0.12.1
     # via matplotlib
-fonttools==4.43.1
+fonttools==4.44.0
     # via matplotlib
 healpy==1.16.6
     # via icecube-skyreader (setup.py)
 idna==3.4
     # via requests
 kiwisolver==1.4.5
     # via matplotlib
-matplotlib==3.8.0
+matplotlib==3.8.1
     # via
     #   healpy
     #   icecube-skyreader (setup.py)
@@ -44,7 +44,7 @@ packaging==23.2
     # via
     #   astropy
     #   matplotlib
-pandas==2.1.1
+pandas==2.1.2
     # via icecube-skyreader (setup.py)
 pillow==10.1.0
     # via matplotlib
@@ -72,5 +72,5 @@ tzdata==2023.3
     # via pandas
 urllib3==2.0.7
     # via requests
-wipac-dev-tools==1.7.0
+wipac-dev-tools==1.7.1
     # via icecube-skyreader (setup.py)

skyreader/plot/contours.py

@@ -0,0 +1,93 @@
+from dataclasses import dataclass
+import healpy  # type: ignore[import]
+import numpy as np
+from typing import ClassVar
+
+
+class CircularContour:
+    @staticmethod
+    def circular_contour(ra_rad: float, dec_rad: float, radius_rad: float, nside: int):
+        """For plotting circular contours on skymaps ra, dec, sigma all
+        expected in radians."""
+        # This transformation is likely reversed later. It can probably be removed.
+        # Likely legacy of old Skymap Scanner design.
+        dec_rad = np.pi / 2.0 - dec_rad
+        step = 1.0 / np.sin(radius_rad) / 10.0
+        # bins * step = 360 deg
+        bins = int(360.0 / step)
+        Theta = np.zeros(bins + 1, dtype=np.double)
+        Phi = np.zeros(bins + 1, dtype=np.double)
+        # Define the contour
+        for j in range(0, bins):
+            # phi runs over 0, 2 pi
+            phi = j * step / 180.0 * np.pi
+            vx = np.cos(phi) * np.sin(ra_rad) * np.sin(radius_rad) + np.cos(ra_rad) * (
+                np.cos(radius_rad) * np.sin(dec_rad)
+                + np.cos(dec_rad) * np.sin(radius_rad) * np.sin(phi)
+            )
+            vy = np.cos(radius_rad) * np.sin(dec_rad) * np.sin(ra_rad) + np.sin(
+                radius_rad
+            ) * (
+                -np.cos(ra_rad) * np.cos(phi)
+                + np.cos(dec_rad) * np.sin(ra_rad) * np.sin(phi)
+            )
+            vz = np.cos(dec_rad) * np.cos(radius_rad) - np.sin(dec_rad) * np.sin(
+                radius_rad
+            ) * np.sin(phi)
+            idx = healpy.vec2pix(nside, vx, vy, vz)
+            Theta[j], Phi[j] = healpy.pix2ang(nside, idx)
+
+        # Close the contour
+        Theta[bins], Phi[bins] = Theta[0], Phi[0]
+
+        return Theta, Phi
+
+
+@dataclass
+class GaussianContour(CircularContour):
+    ra_deg: float
+    dec_deg: float
+    radius_50: float
+    levels: list[float] = [50.0, 90.0]
+    title: str = ""
+    color: str = "m"
+    marker: str = "x"
+    marker_size: int = 20
+
+    CONTOUR_STYLES: ClassVar[dict[float, str]] = {50.0: "-", 90.0: "--"}
+    SCALE_FACTORS: ClassVar[dict[float, float]] = {50.0: 1.177, 90.0: 2.1459}
+
+    @property
+    def dec_rad(self):
+        return np.deg2rad(self.dec_deg)
+
+    @property
+    def ra_rad(self):
+        return np.deg2rad(self.ra_deg)
+
+    @property
+    def sigma_deg(self):
+        return self.radius_50 / self.SCALE_FACTORS[50]
+
+    @property
+    def sigma_rad(self):
+        return np.rad2deg(self.sigma_deg)
+
+    def get_style(self, containment: float):
+        return self.CONTOUR_STYLES[containment]
+
+    def sigma2radius(self, containment: float):
+        # Converts the sigma (standard deviation) to a given containment radius
+        # given a required % of containment.
+        # We use hardcoded values but this can be dynamically calculated with
+        # a chi2 PDF.
+        return self.SCALE_FACTORS[containment]
+
+    def generate_contour(self, nside: int, containment: float):
+        """For plotting circular contours on skymaps ra, dec, sigma all
+        expected in radians."""
+        radius_rad = self.sigma_rad * self.sigma2radius(containment=containment)
+        theta, phi = self.circular_contour(
+            ra_rad=self.ra_rad, dec_rad=self.dec_rad, radius_rad=radius_rad, nside=nside
+        )
+        return theta, phi