Merge pull request #21 from simonsobs/sky_coverage

Sky Coverage Initial Script

Author: Sulla2012

.gitignore

@@ -6,3 +6,4 @@ dist/*
 *.egg-info
 *.fits
 *.hdf
+.coverage

mapcat/core/__init__.py

@@ -0,0 +1,5 @@
+from .core import get_maps_by_coverage
+
+__all__ = [
+    "get_maps_by_coverage",
+]

mapcat/core/core.py

@@ -0,0 +1,52 @@
+from astropy.coordinates import ICRS
+from sqlalchemy import select
+from sqlalchemy.orm import Session
+
+from mapcat.database import DepthOneMapTable
+from mapcat.toolkit.update_sky_coverage import dec_to_index, ra_to_index
+
+
+def get_maps_by_coverage(
+    position: ICRS,
+    session: Session,
+) -> list[DepthOneMapTable]:
+    """
+    Get the depth one maps that cover a given position.
+
+    Parameters
+    ----------
+    position : ICRS
+        The position to query for coverage. Should be in ICRS coordinates.
+    session : Session
+        The database session to use for the query.
+
+    Returns
+    -------
+    session.execute(stmt).scalars().all() : list[DepthOneMapTable]
+        A list of depth one maps that cover the given position.
+
+    Raises
+    ------
+    ValueError
+        If the RA or Dec of the position is out of bounds.
+    """
+    ra = position.ra.deg
+    dec = position.dec.deg
+
+    # These aren't covered since ICRS automatically wraps
+    # values back aground to 0-360 for RA and -90 to 90 for Dec.
+    if ra < 0 or ra > 360:  # pragma: no cover
+        raise ValueError("RA must be between 0 and 360 degrees")
+    if dec < -90 or dec > 90:  # pragma: no cover
+        raise ValueError("Dec must be between -90 and 90 degrees")
+
+    ra_idx = ra_to_index(ra)
+    dec_idx = dec_to_index(dec)
+
+    stmt = (
+        select(DepthOneMapTable)
+        .join(DepthOneMapTable.depth_one_sky_coverage)
+        .filter_by(x=ra_idx, y=dec_idx)
+    )
+
+    return session.execute(stmt).scalars().all()

mapcat/database/atomic_coadd.py

@@ -9,7 +9,7 @@
 if TYPE_CHECKING:
     from .atomic_map import AtomicMapTable
 
-from .links import AtomicMapToCoaddTable, CoaddMapToCoaddTable
+from .links import AtomicMapToCoaddTable, CoaddMapToCoaddTable  # pragma: no cover
 
 
 class AtomicMapCoaddTable(SQLModel, table=True):

mapcat/database/atomic_map.py

@@ -9,7 +9,7 @@
 from .links import AtomicMapToCoaddTable
 
 if TYPE_CHECKING:
-    from .atomic_coadd import AtomicMapCoaddTable
+    from .atomic_coadd import AtomicMapCoaddTable  # pragma: no cover
 
 
 class AtomicMapTable(SQLModel, table=True):

mapcat/database/depth_one_map.py

@@ -6,7 +6,7 @@
 
 from sqlmodel import JSON, Field, Relationship, SQLModel
 
-if TYPE_CHECKING:
+if TYPE_CHECKING:  # pragma: no cover
     from .depth_one_coadd import DepthOneCoaddTable
     from .pipeline_information import PipelineInformationTable
     from .pointing_residual import PointingResidualTable

mapcat/helper.py

@@ -66,7 +66,7 @@ def session(self) -> sessionmaker:
 settings = Settings()
 
 
-def migrate():
+def migrate():  # pragma: no cover
     """
     CLI script to run 'alembic upgrade head' with the correct location.
     """

mapcat/toolkit/plot_tiles.py

@@ -0,0 +1,99 @@
+import argparse as ap
+
+import matplotlib.pyplot as plt
+import numpy as np
+from pixell import enmap
+
+from mapcat.toolkit.update_sky_coverage import *
+
+parser = ap.ArgumentParser()
+parser.add_argument("--imap_path", type=str)
+parser.add_argument("--d1map_path", type=str)
+parser.add_argument("--opath", type=str)
+
+args = parser.parse_args()
+imap_path = args.imap_path
+d1map_path = args.d1map_path
+opath = args.opath
+
+imap = enmap.read_map(str(imap_path))
+
+box = imap.box()
+
+dec_min, ra_max = np.rad2deg(box[0])
+dec_max, ra_min = np.rad2deg(box[1])
+
+pad_low = int((90 + dec_min) * 6 * 2)
+pad_high = int((90 - dec_max) * 6 * 2)
+
+imap = imap[0][::10, ::10]
+
+pad_map = np.pad(
+    imap, ((pad_low, pad_high), (0, 0)), mode="constant", constant_values=0
+)
+del imap
+
+left_limit = pad_map.shape[1]
+right_limit = 0
+top_limit = pad_map.shape[0]
+bottom_limit = 0
+extent = [left_limit, right_limit, bottom_limit, top_limit]
+
+
+plt.imshow(pad_map, vmin=-300, vmax=300, origin="lower", extent=extent)
+plt.vlines(
+    np.arange(0, 360 * 6 * 2, 10 * 6 * 2), ymin=0, ymax=180 * 6 * 2, color="black", lw=1
+)
+plt.hlines(
+    np.arange(0, 180 * 6 * 2, 10 * 6 * 2), xmin=0, xmax=360 * 6 * 2, color="black", lw=1
+)
+plt.xticks(np.arange(0, 360 * 6 * 2, 20 * 6 * 2), labels=np.arange(0, 360, 20))
+plt.yticks(np.arange(0, 180 * 6 * 2, 10 * 6 * 2), labels=np.arange(-90, 90, 10))
+plt.xlabel("RA (degrees)")
+plt.ylabel("Dec (degrees)")
+
+d1map = enmap.read_map(str(d1map_path))
+coverage_tiles = get_sky_coverage(d1map)
+
+d1box = d1map.box()
+
+d1dec_min, d1ra_max = np.rad2deg(d1box[0])
+d1dec_max, d1ra_min = np.rad2deg(d1box[1])
+
+d1pad_low_dec = int((90 + d1dec_min) * 6 * 2)
+d1pad_high_dec = int((90 - d1dec_max) * 6 * 2)
+
+d1pad_low_ra = int((180 + d1ra_min) * 6 * 2)
+d1pad_high_ra = int((180 - d1ra_max) * 6 * 2)
+
+d1map = d1map[0][::10, ::10]
+d1pad_map = np.pad(
+    d1map,
+    ((d1pad_low_dec, d1pad_high_dec), (d1pad_high_ra, d1pad_low_ra)),
+    mode="constant",
+    constant_values=0,
+)
+
+plt.imshow(
+    d1pad_map,
+    vmin=-300,
+    vmax=300,
+    origin="lower",
+    alpha=0.5,
+    cmap="seismic",
+    extent=extent,
+)
+
+for tile in coverage_tiles:
+    plt.gca().add_patch(
+        plt.Rectangle(
+            (tile[0] * 10 * 6 * 2, tile[1] * 10 * 6 * 2),
+            10 * 6 * 2,
+            10 * 6 * 2,
+            fill=False,
+            edgecolor="red",
+            lw=2,
+        )
+    )
+
+plt.savefig(opath + "act_coverage.png", dpi=300)

mapcat/toolkit/reset.py

@@ -123,8 +123,7 @@ def main():
         type=float,
         default=None,
         help=(
-            "Only reset entries for maps whose ctime is >= START_TIME "
-            "(unix timestamp)."
+            "Only reset entries for maps whose ctime is >= START_TIME (unix timestamp)."
         ),
     )
 
@@ -133,8 +132,7 @@ def main():
         type=float,
         default=None,
         help=(
-            "Only reset entries for maps whose ctime is <= END_TIME "
-            "(unix timestamp)."
+            "Only reset entries for maps whose ctime is <= END_TIME (unix timestamp)."
         ),
     )
 

mapcat/toolkit/update_sky_coverage.py

@@ -25,6 +25,88 @@ def resolve_tmap(d1table: DepthOneMapTable) -> Path:
     return settings.depth_one_parent / d1table.mean_time_path
 
 
+def index_to_skybox(ra_idx: int, dec_idx: int) -> np.ndarray:
+    """
+    Convert a sky coverage tile index to a sky box in radians
+
+    Parameters
+    ----------
+    ra_idx : int
+        The RA index of the sky coverage tile
+    dec_idx : int
+        The Dec index of the sky coverage tile
+
+    Returns
+    -------
+    skybox : np.ndarray
+        A 2x2 array containing the corners of the sky box in radians, in the format [[dec_min, ra_max], [dec_max, ra_min]]
+    """
+    ra_min = ra_idx * 10
+    ra_max = ra_min + 10
+    dec_min = (dec_idx - 9) * 10
+    dec_max = dec_min + 10
+
+    return np.array(
+        [
+            [np.deg2rad(dec_min), np.deg2rad(ra_max)],
+            [np.deg2rad(dec_max), np.deg2rad(ra_min)],
+        ]
+    )
+
+
+def ra_to_index(ra: float) -> int:
+    """
+    Convert an ra in degrees to a sky coverage tile index
+
+    Parameters
+    ----------
+    ra : float
+        The ra in degrees to convert
+
+    Returns
+    -------
+    idx : int
+        The sky coverage tile index corresponding to the input ra
+    """
+    return int(np.floor(ra / 10))
+
+
+def dec_to_index(dec: float) -> int:
+    """
+    Convert a dec in degrees to a sky coverage tile index
+
+    Parameters
+    ----------
+    dec : float
+        The dec in degrees to convert
+
+    Returns
+    -------
+    idx : int
+        The sky coverage tile index corresponding to the input dec
+    """
+    return int(np.floor(dec / 10)) + 9
+
+
+def _ra_to_index_pixell(ra: float) -> int:
+    """
+    Convert an ra in degrees to a sky coverage tile index using the
+    pixell convention where -180 < ra < 180. You should probably
+    not ever touch this function.
+
+    Parameters
+    ----------
+    ra : float
+        The ra in degrees to convert
+
+    Returns
+    -------
+    idx : int
+        The sky coverage tile index corresponding to the input ra
+    """
+    return int(np.floor(ra / 10)) + 18
+
+
 def get_sky_coverage(tmap: enmap.ndmap) -> list:
     """
     Given the time map of a depth1 map, return the list
@@ -49,27 +131,27 @@ def get_sky_coverage(tmap: enmap.ndmap) -> list:
     dec_max = np.ceil(dec_max / 10) * 10
     ra_min = np.floor(ra_min / 10) * 10
     ra_max = np.ceil(ra_max / 10) * 10
+    ra_min += 180  # Convert from pixel standard to normal RA convention
+    ra_max += 180
 
     ras = np.arange(ra_min, ra_max, 10)
     decs = np.arange(dec_min, dec_max, 10)
 
     ra_idx = []
-    dec_id = []
+    dec_idx = []
 
     for ra in ras:
         for dec in decs:
-            skybox = np.array(
-                [
-                    [np.deg2rad(dec), np.deg2rad(ra + 10)],
-                    [np.deg2rad(dec + 10), np.deg2rad(ra)],
-                ]
-            )
+            ra_id = ra_to_index(ra)
+            dec_id = dec_to_index(dec)
+            skybox = index_to_skybox(ra_id, dec_id)
+            skybox[..., 1] -= np.pi  # Convert from standard RA to pixell convention
             submap = enmap.submap(tmap, skybox)
             if np.any(submap):
-                ra_idx.append(int(ra / 10))
-                dec_id.append(int(dec / 10) + 9)
+                ra_idx.append(ra_id)
+                dec_idx.append(dec_id)
 
-    return list(zip(ra_idx, dec_id))
+    return list(zip(ra_idx, dec_idx))
 
 
 def coverage_from_depthone(d1table: DepthOneMapTable) -> list[SkyCoverageTable]: