diff --git a/.gitignore b/.gitignore index 54ee18a..db407ce 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ dist/* *.egg-info *.fits *.hdf +.coverage diff --git a/mapcat/core/__init__.py b/mapcat/core/__init__.py new file mode 100644 index 0000000..3bee488 --- /dev/null +++ b/mapcat/core/__init__.py @@ -0,0 +1,5 @@ +from .core import get_maps_by_coverage + +__all__ = [ + "get_maps_by_coverage", +] diff --git a/mapcat/core/core.py b/mapcat/core/core.py new file mode 100644 index 0000000..e324388 --- /dev/null +++ b/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() diff --git a/mapcat/database/atomic_coadd.py b/mapcat/database/atomic_coadd.py index fad6a61..73653e9 100644 --- a/mapcat/database/atomic_coadd.py +++ b/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): diff --git a/mapcat/database/atomic_map.py b/mapcat/database/atomic_map.py index 57304a7..08753a3 100644 --- a/mapcat/database/atomic_map.py +++ b/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): diff --git a/mapcat/database/depth_one_map.py b/mapcat/database/depth_one_map.py index aa32033..bb8f7f6 100644 --- a/mapcat/database/depth_one_map.py +++ b/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 diff --git a/mapcat/helper.py b/mapcat/helper.py index 6efc7d6..29d348a 100644 --- a/mapcat/helper.py +++ b/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. """ diff --git a/mapcat/toolkit/plot_tiles.py b/mapcat/toolkit/plot_tiles.py new file mode 100644 index 0000000..71467e2 --- /dev/null +++ b/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) diff --git a/mapcat/toolkit/reset.py b/mapcat/toolkit/reset.py index 8cdcc6d..d055bd5 100644 --- a/mapcat/toolkit/reset.py +++ b/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)." ), ) diff --git a/mapcat/toolkit/update_sky_coverage.py b/mapcat/toolkit/update_sky_coverage.py index ec77a19..3679ad0 100644 --- a/mapcat/toolkit/update_sky_coverage.py +++ b/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]: diff --git a/pyproject.toml b/pyproject.toml index e09f6ad..118d222 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,7 @@ dependencies = [ [project.scripts] actingest = "mapcat.toolkit.act:main" mapcatmigrate = "mapcat.helper:migrate" +updatesky = "mapcat.toolkit.update_sky_coverage:main" mapcatreset = "mapcat.toolkit.reset:main" [project.optional-dependencies] diff --git a/tests/test_act.py b/tests/test_act.py index e16cbb0..a89cba5 100644 --- a/tests/test_act.py +++ b/tests/test_act.py @@ -2,11 +2,16 @@ import os from pathlib import Path +import astropy.units as u +import numpy as np import pytest import requests +from astropy.coordinates import ICRS +from pixell import enmap from sqlalchemy import create_engine from sqlalchemy.orm import sessionmaker +from mapcat.core import get_maps_by_coverage from mapcat.database import DepthOneMapTable from mapcat.toolkit import act, update_sky_coverage @@ -27,86 +32,86 @@ cov_mapping = { "1505603190.0": [ - (-4, 3), - (-4, 4), - (-4, 6), - (-4, 7), - (-3, 3), - (-3, 4), - (-3, 5), - (-3, 6), - (-3, 7), - (-2, 3), - (-2, 4), - (-2, 5), - (-2, 6), - (-2, 7), - (-1, 3), - (-1, 4), - (-1, 5), - (-1, 6), - (-1, 7), - (0, 3), - (0, 4), - (0, 5), - (0, 6), - (0, 7), - (1, 3), - (1, 4), - (1, 5), - (1, 6), - (1, 7), - (2, 3), - (2, 4), - (2, 5), - (2, 6), - (2, 7), - (3, 3), - (3, 4), - (3, 5), - (3, 6), - (3, 7), - (4, 3), - (4, 4), - (4, 5), - (4, 6), - (4, 7), - (5, 3), - (5, 4), - (5, 5), - (5, 6), - (5, 7), - (6, 3), - (6, 4), - (6, 5), - (6, 6), - (6, 7), - (7, 3), - (7, 4), - (7, 5), - (7, 6), - (7, 7), - (8, 3), - (8, 4), - (8, 5), - (8, 6), - (8, 7), - (9, 3), - (9, 4), - (9, 5), - (9, 6), - (9, 7), - (10, 3), - (10, 4), - (10, 5), - (10, 6), - (10, 7), - (11, 3), - (11, 4), - (11, 5), - (11, 6), + (14, 3), + (14, 4), + (14, 6), + (14, 7), + (15, 3), + (15, 4), + (15, 5), + (15, 6), + (15, 7), + (16, 3), + (16, 4), + (16, 5), + (16, 6), + (16, 7), + (17, 3), + (17, 4), + (17, 5), + (17, 6), + (17, 7), + (18, 3), + (18, 4), + (18, 5), + (18, 6), + (18, 7), + (19, 3), + (19, 4), + (19, 5), + (19, 6), + (19, 7), + (20, 3), + (20, 4), + (20, 5), + (20, 6), + (20, 7), + (21, 3), + (21, 4), + (21, 5), + (21, 6), + (21, 7), + (22, 3), + (22, 4), + (22, 5), + (22, 6), + (22, 7), + (23, 3), + (23, 4), + (23, 5), + (23, 6), + (23, 7), + (24, 3), + (24, 4), + (24, 5), + (24, 6), + (24, 7), + (25, 3), + (25, 4), + (25, 5), + (25, 6), + (25, 7), + (26, 3), + (26, 4), + (26, 5), + (26, 6), + (26, 7), + (27, 3), + (27, 4), + (27, 5), + (27, 6), + (27, 7), + (28, 3), + (28, 4), + (28, 5), + (28, 6), + (28, 7), + (29, 3), + (29, 4), + (29, 5), + (29, 6), ], - "1505646390.0": [(4, 3), (4, 4), (4, 5), (5, 3), (5, 4), (5, 5)], + "1505646390.0": [(22, 3), (22, 4), (22, 5), (23, 3), (23, 4), (23, 5)], } @@ -222,6 +227,7 @@ def test_act(database_sessionmaker, downloaded_data_file): def test_sky_coverage(database_sessionmaker, downloaded_data_file): + # I'm not sure we need this test any more, test_sky_coverage_2 is better args = ap.Namespace( glob="*/*_map.fits", relative_to=downloaded_data_file, @@ -243,3 +249,59 @@ def test_sky_coverage(database_sessionmaker, downloaded_data_file): ), # These should be ints, idk why I have to cast them (from str) int(cov.y), ) in cov_mapping[str(d1map.ctime)] + + with database_sessionmaker() as session: + maps = session.query(DepthOneMapTable).all() + for map in maps: + # Clean up + session.delete(map) + + session.commit() + + +def test_sky_coverage_2(database_sessionmaker, downloaded_data_file): + args = ap.Namespace( + glob="*/*_map.fits", + relative_to=downloaded_data_file, + telescope="act", + ) + + act.core(session=database_sessionmaker, args=args) + + update_sky_coverage.core(session=database_sessionmaker) + + d1maps = act.glob(args.glob, args.relative_to, args.telescope) + with database_sessionmaker() as session: + for d1map in d1maps: + cur_map = enmap.read_map(str(downloaded_data_file) + "/" + d1map.map_path) + nonzero_radec = cur_map.pix2sky(np.where(cur_map[0] != 0)) + idx = np.linspace(0, len(nonzero_radec) - 1, 1000) + idx = np.round(idx).astype(int) + nonzero_radec = nonzero_radec.T[ + idx + ] # Only test a subset of the nonzero pixels to speed up the test + for pix in nonzero_radec: + coord = ICRS( + (pix[1] + np.pi) * u.rad, pix[0] * u.rad + ) # Convert from pixel standard to normal RA convention + return_d1map = get_maps_by_coverage(coord, session) + assert d1map.map_name in [m.map_name for m in return_d1map] + + ra, dec = 0, 0 # Test a point outside the coverage, should return an empty list + coord = ICRS(ra * u.rad, dec * u.rad) + return_d1map = get_maps_by_coverage(coord, session) + assert len(return_d1map) == 0 + + with database_sessionmaker() as session: + maps = session.query(DepthOneMapTable).all() + for map in maps: + # Clean up + session.delete(map) + + session.commit() + + +def test_ra_to_index_pixell(): + assert update_sky_coverage._ra_to_index_pixell(-180) == 0 + assert update_sky_coverage._ra_to_index_pixell(170) == 35 + assert update_sky_coverage._ra_to_index_pixell(0) == 18 diff --git a/tests/test_reset.py b/tests/test_reset.py index a2fbcb1..2650631 100644 --- a/tests/test_reset.py +++ b/tests/test_reset.py @@ -19,9 +19,7 @@ def run_migration(database_path: str): from alembic import command from alembic.config import Config - alembic_cfg = Config( - str(Path(__file__).parent.parent / "mapcat" / "alembic.ini") - ) + alembic_cfg = Config(str(Path(__file__).parent.parent / "mapcat" / "alembic.ini")) database_url = f"sqlite:///{database_path}" alembic_cfg.set_main_option("sqlalchemy.url", database_url) command.upgrade(alembic_cfg, "head")