diff --git a/README.md b/README.md index c41e05a..fd5240b 100644 --- a/README.md +++ b/README.md @@ -58,6 +58,23 @@ actingest --relative-to=/path/to/maps --glob=*/*_map.fits --telescope=act ``` More information on the parameters is available through `actingest -h`. +### Sky coverage + +After ingestion, generate the 10-degree coverage tiles with `updatesky`. +Use the same `--relative-to` directory as `actingest`, or configure +`MAPCAT_DEPTH_ONE_PARENT` as the map root. A tile is +covered when it contains a finite, nonzero pixel center. RA wraps into +0–360 degrees; tile indices are 0–35 in RA and 0–17 in declination. +Signed RA in ACT FITS files represents the same celestial longitude as +unsigned RA. + +Coverage can be regenerated using the following command: +``` +updatesky --replace +``` +This recomputes all catalog maps and replaces their coverage in one +transaction. Without `--replace`, maps with existing coverage are skipped. + Registering new Maps -------------------- diff --git a/mapcat/database/sky_coverage.py b/mapcat/database/sky_coverage.py index 3906f6f..a05dc67 100644 --- a/mapcat/database/sky_coverage.py +++ b/mapcat/database/sky_coverage.py @@ -13,7 +13,7 @@ class SkyCoverageTable(SQLModel, table=True): """ Table for tracking sky coverage patches with non-zero overlap with a given depth one map. - x and y are 0->36 and 0-18 respectively for CAR patches with 10x10 degrees each. + x and y are 0..35 and 0..17 respectively for CAR patches with 10x10 degrees each. Attributes ---------- diff --git a/mapcat/toolkit/plot_tiles.py b/mapcat/toolkit/plot_tiles.py index dab8d44..b8343a0 100644 --- a/mapcat/toolkit/plot_tiles.py +++ b/mapcat/toolkit/plot_tiles.py @@ -1,99 +1,104 @@ import argparse as ap +from itertools import pairwise +from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pixell import enmap -from mapcat.toolkit.update_sky_coverage import get_sky_coverage - -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, convention="ACT") - -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, -) +from mapcat.toolkit.update_sky_coverage import _separable_car_axes, get_sky_coverage + + +def _show_map(map_data, step=10, **kwargs): + """Display sampled CAR pixels in celestial coordinates, splitting at RA=0. + + Raises + ------ + ValueError + If the map is not an unrotated equatorial CAR map. + """ + axes = _separable_car_axes(map_data) + if axes is None: + raise ValueError("plot_tiles requires an unrotated equatorial CAR map") + dec, ra = axes + # Use the original pixel centers: striding an enmap changes its WCS centers. + dec, ra = dec[::step], ra[::step] % 360 + pixels = np.asarray(map_data.preflat[0])[::step, ::step] + matrix = map_data.wcs.wcs.get_pc() + dra, ddec = map_data.wcs.wcs.cdelt * np.diag(matrix) * step + boundaries = np.r_[0, np.flatnonzero(np.abs(np.diff(ra)) > 180) + 1, len(ra)] + for start, stop in pairwise(boundaries): + extent = [ + ra[start] - dra / 2, + ra[stop - 1] + dra / 2, + dec[0] - ddec / 2, + dec[-1] + ddec / 2, + ] + plt.imshow(pixels[:, start:stop], origin="lower", extent=extent, **kwargs) + # A pixel footprint straddling RA=0 also appears at the other sky edge. + for shift in ([360] if min(extent[:2]) < 0 else []) + ( + [-360] if max(extent[:2]) > 360 else [] + ): + plt.imshow( + pixels[:, start:stop], + origin="lower", + extent=[extent[0] + shift, extent[1] + shift, *extent[2:]], + **kwargs, + ) + + +def main(): + parser = ap.ArgumentParser() + parser.add_argument("--imap_path", type=str, required=True) + parser.add_argument("--d1map_path", type=str, required=True) + parser.add_argument("--opath", type=str, required=True) + args = parser.parse_args() + + # Preserve the existing palette and optional notebook style. + try: + import socolors # noqa: F401 + except ImportError: + pass + if "notebook" in plt.style.available: + plt.style.use("notebook") + + imap = enmap.read_map(str(args.imap_path), preflat=True, sel=0) + _show_map(imap, vmin=-200, vmax=200, zorder=-1000, cmap="grey") + del imap + + plt.vlines(np.arange(0, 360, 10), ymin=-90, ymax=90, color="black", lw=1) + plt.hlines(np.arange(-90, 90, 10), xmin=0, xmax=360, color="black", lw=1) + plt.xticks(np.arange(0, 360, 20), labels=np.arange(0, 360, 20)) + plt.yticks(np.arange(-90, 90, 10), labels=np.arange(-90, 90, 10)) + plt.xlabel("RA (degrees)") + plt.ylabel("Dec (degrees)") + + d1map = enmap.read_map(str(args.d1map_path)) + coverage_tiles = get_sky_coverage(d1map) + sampled = np.asarray(d1map.preflat[0])[::10, ::10] + observed = sampled[np.isfinite(sampled) & (sampled != 0)] + vmin, vmax = np.percentile(observed, (1, 99)) if observed.size else (0, 1) + _show_map(d1map, vmin=vmin, vmax=vmax, alpha=0.9, cmap="twilight_shifted") + + for tile in coverage_tiles: + plt.gca().add_patch( + plt.Rectangle( + (tile[0] * 10, tile[1] * 10 - 90), + 10, + 10, + fill=False, + edgecolor="C0", + lw=2, + ) + ) -plt.imshow( - d1pad_map, - vmin=-300, - vmax=300, - origin="lower", - alpha=0.5, - cmap="seismic", - extent=extent, -) + plt.xlim(360, 0) + plt.ylim(-90, 90) + opath = Path(args.opath) + opath.mkdir(parents=True, exist_ok=True) + plt.savefig(opath / "act_coverage.png", dpi=300) + plt.close() -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) +if __name__ == "__main__": + main() diff --git a/mapcat/toolkit/update_sky_coverage.py b/mapcat/toolkit/update_sky_coverage.py index b2d1139..494fdc8 100644 --- a/mapcat/toolkit/update_sky_coverage.py +++ b/mapcat/toolkit/update_sky_coverage.py @@ -1,14 +1,17 @@ +import argparse +from collections.abc import Iterator from pathlib import Path import numpy as np from pixell import enmap +from tqdm import tqdm from mapcat.database.depth_one_map import DepthOneMapTable from mapcat.database.sky_coverage import SkyCoverageTable from mapcat.helper import settings -def resolve_tmap(d1table: DepthOneMapTable) -> Path: +def resolve_tmap(d1table: DepthOneMapTable, *, relative_to: Path | None = None) -> Path: """ Resolve the local path to a tmap from a d1 table. @@ -16,13 +19,24 @@ def resolve_tmap(d1table: DepthOneMapTable) -> Path: ---------- d1table : DepthOneMapTable The depth one map table to resolve the tmap for + relative_to : Path, optional + Base directory used when ingesting the map. Defaults to the configured + depth_one_parent. Returns ------- - ettings.depth_one_parent / d1table.mean_time_path : Path + path : Path The local path to the tmap for the depth one map + + Raises + ------ + ValueError + If the map has no mean time path. """ - return settings.depth_one_parent / d1table.mean_time_path + if d1table.mean_time_path is None: + raise ValueError(f"No mean time map available for {d1table.map_name}") + parent = settings.depth_one_parent if relative_to is None else Path(relative_to) + return parent / d1table.mean_time_path def index_to_skybox(ra_idx: int, dec_idx: int) -> np.ndarray: @@ -67,8 +81,15 @@ def ra_to_index(ra: float) -> int: ------- idx : int The sky coverage tile index corresponding to the input ra + + Raises + ------ + ValueError + If RA is not finite. """ - return int(np.floor(ra / 10)) + if not np.isfinite(ra): + raise ValueError("RA must be finite") + return int(_tile_floor((ra % 360) / 10)) % 36 def dec_to_index(dec: float) -> int: @@ -84,111 +105,185 @@ def dec_to_index(dec: float) -> int: ------- idx : int The sky coverage tile index corresponding to the input dec + + Raises + ------ + ValueError + If Dec is nonfinite or outside [-90, 90]. """ - return int(np.floor(dec / 10)) + 9 + if not np.isfinite(dec) or not -90 <= dec <= 90: + raise ValueError("Dec must be between -90 and 90 degrees") + return min(int(_tile_floor((dec + 90) / 10)), 17) + + +def _tile_floor(values): + """Stabilize exact tile boundaries against FITS WCS floating-point noise.""" + nearest = np.rint(values) + return np.floor(np.where(np.abs(values - nearest) < 1e-10, nearest, values)) + +def _observed_pixel_coordinates( + tmap: enmap.ndmap, +) -> Iterator[tuple[np.ndarray, np.ndarray]]: + """Yield observed pixel rows and columns, processing 128 rows at a time. -def _ra_to_index_pixell(ra: float) -> int: + A pixel is observed if any component has a finite, nonzero value. + Returned row numbers refer to the full map, not the current block. """ - 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. + for first_row in range(0, tmap.shape[-2], 128): + block = np.asarray(tmap[..., first_row : first_row + 128, :]) + observed = np.isfinite(block) & (block != 0) + component_axes = tuple(range(observed.ndim - 2)) + if component_axes: + observed = np.any(observed, axis=component_axes) - Parameters - ---------- - ra : float - The ra in degrees to convert + rows, columns = np.nonzero(observed) + if rows.size: + yield rows + first_row, columns - Returns - ------- - idx : int - The sky coverage tile index corresponding to the input ra + +def _separable_car_axes( + tmap: enmap.ndmap, +) -> tuple[np.ndarray, np.ndarray] | None: + """Cache Dec per row and RA per column for unrotated equatorial CAR maps. + + Each sky tile is a rectangle of rows and columns for these maps, even + when either pixel axis is reversed. Return None for other geometries. """ - return int(np.floor(ra / 10)) + 18 + wcs = tmap.wcs.wcs + if list(wcs.ctype) != ["RA---CAR", "DEC--CAR"] or wcs.crval[1] != 0: + return None + + # Off-diagonal matrix terms mix rows and columns, preventing this shortcut. + pixel_matrix = wcs.get_pc() + if not np.all(pixel_matrix == np.diag(np.diag(pixel_matrix))): + return None + + nrows, ncolumns = tmap.shape[-2:] + dec_by_row = np.rad2deg( + enmap.pix2sky( + tmap.shape, tmap.wcs, [np.arange(nrows), np.zeros(nrows)], safe=False + )[0] + ) + ra_by_column = np.rad2deg( + enmap.pix2sky( + tmap.shape, tmap.wcs, [np.zeros(ncolumns), np.arange(ncolumns)], safe=False + )[1] + ) + return dec_by_row, ra_by_column -def get_sky_coverage(tmap: enmap.ndmap, convention: str = "standard") -> list: +def _sky_coordinates_to_tiles( + dec_degrees: np.ndarray, + ra_degrees: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Convert valid celestial coordinates into RA and Dec tile indices. + + RA wraps around the sky. Dec is offset by 90 degrees to index from the + south pole; clipping keeps the north pole in the final tile. Boundary + rounding uses the same tolerance as the scalar coordinate helpers. """ - Given the time map of a depth1 map, return the list - of sky coverage tiles that cover that map + valid = ( + np.isfinite(ra_degrees) & np.isfinite(dec_degrees) & (np.abs(dec_degrees) <= 90) + ) + ra_tiles = _tile_floor((ra_degrees[valid] % 360) / 10).astype(int) % 36 + dec_tiles = np.clip(_tile_floor((dec_degrees[valid] + 90) / 10).astype(int), 0, 17) + return ra_tiles, dec_tiles + + +def _coverage_for_nonseparable_wcs(tmap: enmap.ndmap) -> list[tuple[int, int]]: + """Handle projections where a sky tile is not a pixel-axis rectangle.""" + covered_tiles = np.zeros((36, 18), dtype=bool) + for rows, columns in _observed_pixel_coordinates(tmap): + dec_degrees, ra_degrees = np.rad2deg( + enmap.pix2sky(tmap.shape, tmap.wcs, [rows, columns], safe=False) + ) + ra_tiles, dec_tiles = _sky_coordinates_to_tiles(dec_degrees, ra_degrees) + covered_tiles[ra_tiles, dec_tiles] = True + + ra_tiles, dec_tiles = np.nonzero(covered_tiles) + return [(int(ra), int(dec)) for ra, dec in zip(ra_tiles, dec_tiles)] + + +def get_sky_coverage(tmap: enmap.ndmap) -> list[tuple[int, int]]: + """Return sorted 10-degree tiles containing finite, nonzero pixel centers. + + For unrotated equatorial CAR maps, including ACT, loop over RA and Dec + tiles and check their pixel rectangles. Select rows and columns by their + WCS coordinates rather than submap rounding or assumed axis directions. + Leading component axes are combined: any observed component marks a pixel. + Other projections use a pixel-coordinate fallback. Parameters ---------- - tmap : enmap.enmap - The time map of the depth-one map. Pixels that were observed have non-zero values. - convention : str, optional - The coordinate convention to use. Default is "standard", corresponding to 0 < ra < 360. - If ACT is specified, the convention is -180 < ra < 180. + tmap : enmap.ndmap + Coverage map with celestial WCS and at least two spatial axes. Returns ------- - tiles : list - A list of sky coverage tiles that cover the map + list + Unique (RA, Dec) tile indices, with 0 <= RA < 36 and 0 <= Dec < 18. Raises - -------- + ------ ValueError - If the convention is not 'standard' or 'ACT'. + If the map has fewer than two axes. """ - if convention not in ["standard", "ACT"]: - raise ValueError("Invalid convention. Must be 'standard' or 'ACT'.") - box = tmap.box() - - dec_min, ra_max = np.rad2deg(box[0]) - dec_max, ra_min = np.rad2deg(box[1]) - - dec_min = np.floor(dec_min / 10) * 10 - dec_max = np.ceil(dec_max / 10) * 10 - ra_min = np.floor(ra_min / 10) * 10 - ra_max = np.ceil(ra_max / 10) * 10 - if convention == "ACT": - 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_idx = [] + if tmap.ndim < 2: + raise ValueError("Coverage maps must have at least two spatial axes") + + coordinate_axes = _separable_car_axes(tmap) + if coordinate_axes is None: + return _coverage_for_nonseparable_wcs(tmap) + dec_by_row, ra_by_column = coordinate_axes + + # Assign pixel centers to tiles. Modulo wraps RA without rotating the sky. + valid_ra = np.isfinite(ra_by_column) + valid_dec = np.isfinite(dec_by_row) & (np.abs(dec_by_row) <= 90) + column_tiles = np.full(ra_by_column.shape, -1, dtype=int) + row_tiles = np.full(dec_by_row.shape, -1, dtype=int) + column_tiles[valid_ra] = ( + _tile_floor((ra_by_column[valid_ra] % 360) / 10).astype(int) % 36 + ) + row_tiles[valid_dec] = np.clip( + _tile_floor((dec_by_row[valid_dec] + 90) / 10).astype(int), 0, 17 + ) - for ra in ras: - for dec in decs: - ra_id = ra_to_index(ra) - dec_id = dec_to_index(dec) - skybox = index_to_skybox(ra_id, dec_id) - if convention == "ACT": - skybox[..., 1] -= np.pi # Convert from standard RA to pixell convention - submap = enmap.submap(tmap, skybox) - if np.any(submap): - ra_idx.append(ra_id) - dec_idx.append(dec_id) + tiles = [] + for ra_idx in np.unique(column_tiles[valid_ra]): + columns = np.flatnonzero(column_tiles == ra_idx) + for dec_idx in np.unique(row_tiles[valid_dec]): + rows = np.flatnonzero(row_tiles == dec_idx) + submap = np.asarray(tmap[..., rows[:, None], columns]) + if np.any(np.isfinite(submap) & (submap != 0)): + tiles.append((int(ra_idx), int(dec_idx))) - return list(zip(ra_idx, dec_idx)) + return tiles def coverage_from_depthone( - d1table: DepthOneMapTable, convention: str = "standard" + d1table: DepthOneMapTable, *, relative_to: Path | None = None ) -> list[SkyCoverageTable]: """ Get the list of sky coverage tiles that cover a given depth one map Parameters ---------- - d1map : DepthOneMapTable + d1table : DepthOneMapTable The depth one map to get the sky coverage for - convention : str, optional - The coordinate convention to use. Default is "standard", corresponding to 0 < ra < 360. - If ACT is specified, the convention is -180 < ra < 180. + relative_to : Path, optional + Base directory used when ingesting the map. Defaults to the configured + depth_one_parent. Returns ------- tiles : list[SkyCoverageTable] A list of sky coverage tiles that cover the map """ - tmap_path = resolve_tmap(d1table) + tmap_path = resolve_tmap(d1table, relative_to=relative_to) tmap = enmap.read_map(str(tmap_path)) - coverage_tiles = get_sky_coverage(tmap, convention=convention) + coverage_tiles = get_sky_coverage(tmap) return [ SkyCoverageTable(x=tile[0], y=tile[1], map_id=d1table.map_id) @@ -196,7 +291,8 @@ def coverage_from_depthone( ] -def core(session, convention: str = "standard"): + +def core(session, *, replace: bool = False, relative_to: Path | None = None, progress: bool = False): """ Core function for updating the sky coverage table. For each depth one map that does not have any associated sky coverage tiles, compute the sky coverage tiles and add them to the database. @@ -204,25 +300,43 @@ def core(session, convention: str = "standard"): ---------- session : sessionmaker A SQLAlchemy sessionmaker to use for database access. - convention : str, optional - The coordinate convention to use. Default is "standard", corresponding to 0 < ra < 360. - If ACT is specified, the convention is -180 < ra < 180. + replace : bool, optional + Recompute existing coverage as well, replacing it in one transaction. + relative_to : Path, optional + Base directory used when ingesting the maps. Defaults to the configured + depth_one_parent. + progress : bool, optional + Show a progress bar for processing maps. """ with session() as cur_session: - d1maps = ( - cur_session.query(DepthOneMapTable) - .outerjoin( + query = cur_session.query(DepthOneMapTable) + if not replace: + query = query.outerjoin( SkyCoverageTable, SkyCoverageTable.map_id == DepthOneMapTable.map_id - ) - .filter(SkyCoverageTable.map_id.is_(None)) - .all() - ) - for d1map in d1maps: - SkyCov = coverage_from_depthone(d1map, convention=convention) + ).filter(SkyCoverageTable.map_id.is_(None)) + for d1map in tqdm(query.all(), disable=not progress, desc="Updating sky coverage"): + SkyCov = coverage_from_depthone(d1map, relative_to=relative_to) + if replace: + cur_session.query(SkyCoverageTable).filter_by( + map_id=d1map.map_id + ).delete(synchronize_session=False) cur_session.add_all(SkyCov) cur_session.commit() def main(): - core(session=settings.session) + parser = argparse.ArgumentParser( + description="Generate sky coverage from time maps." + ) + parser.add_argument( + "--replace", action="store_true", help="Recompute and replace existing coverage" + ) + parser.add_argument( + "--relative-to", type=Path, help="Base directory used for ACT ingestion" + ) + parser.add_argument( + "--progress", action="store_true", help="Show progress bar for processing maps" + ) + args = parser.parse_args() + core(session=settings.session, replace=args.replace, relative_to=args.relative_to, progress=args.progress) diff --git a/tests/test_act.py b/tests/test_act.py index 60a2741..4ddc660 100644 --- a/tests/test_act.py +++ b/tests/test_act.py @@ -16,19 +16,17 @@ from mapcat.database import DepthOneMapTable from mapcat.toolkit import act, update_sky_coverage +BASE_URL = "https://object-arbutus.alliancecan.ca/f620008d8888477e9fc9e5dca514fbc3:sodacan-public/maps/act/15056" +INFO_BASE_URL = "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056" + +# The public bucket supplies FITS maps, but not the ingestion metadata HDFs. DATA_URLS = [ - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505603190_pa4_f150_info.hdf", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505603190_pa4_f150_ivar.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505603190_pa4_f150_kappa.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505603190_pa4_f150_map.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505603190_pa4_f150_rho.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505603190_pa4_f150_time.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505646390_pa6_f150_info.hdf", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505646390_pa6_f150_ivar.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505646390_pa6_f150_kappa.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505646390_pa6_f150_map.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505646390_pa6_f150_rho.fits", - "https://g-0a470a.6b7bd8.0ec8.data.globus.org/act_dr6/dr6.02/depth1/depth1_maps/15056/depth1_1505646390_pa6_f150_time.fits", + f"{INFO_BASE_URL}/depth1_1505603190_pa4_f150_info.hdf", + f"{BASE_URL}/depth1_1505603190_pa4_f150_map.fits", + f"{BASE_URL}/depth1_1505603190_pa4_f150_time.fits", + f"{INFO_BASE_URL}/depth1_1505646390_pa6_f150_info.hdf", + f"{BASE_URL}/depth1_1505646390_pa6_f150_map.fits", + f"{BASE_URL}/depth1_1505646390_pa6_f150_time.fits", ] cov_mapping = { @@ -115,6 +113,12 @@ "1505646390.0": [(22, 3), (22, 4), (22, 5), (23, 3), (23, 4), (23, 5)], } +# The historical expectations used RA + 180 degrees; use celestial RA. +cov_mapping = { + ctime: sorted(((x + 18) % 36, y) for x, y in tiles) + for ctime, tiles in cov_mapping.items() +} + def run_migration(database_path: str): """ @@ -131,7 +135,7 @@ def run_migration(database_path: str): command.upgrade(alembic_cfg, "head") -@pytest.fixture(scope="session", autouse=True) +@pytest.fixture(autouse=True) def database_sessionmaker(tmp_path_factory): """ Create a temporary SQLite database for testing. @@ -150,12 +154,14 @@ def database_sessionmaker(tmp_path_factory): yield sessionmaker(bind=engine, expire_on_commit=False) + engine.dispose() + # Clean up the database (don't do this in case we want to inspect) database_path.unlink() @pytest.fixture(scope="session") -def downloaded_data_file(request): +def downloaded_data_file(request, pytestconfig): """ Fixture to download depth 1 maps for testing. @@ -169,37 +175,39 @@ def downloaded_data_file(request): Path to the downloaded file """ + cache_dir = Path(pytestconfig.cache._cachedir) / "d1maps" for url in DATA_URLS: - cache_dir = Path("../.pytest_cache/d1maps") filename = os.path.basename(url) subdir = os.path.basename(os.path.dirname(url)) file_path = cache_dir / subdir / filename - # Check to see if each file is already downloaded in the cache, if so skip downloading it again - cached_file = request.config.cache.get( - "downloaded_file_" + os.path.basename(url), None - ) - if cached_file and os.path.exists(cached_file): + # The directory is the source of truth; cached strings may point to + # a previous working directory after moving or restoring the cache. + if file_path.is_file(): continue # Download the file - response = requests.get(url) - response.raise_for_status() - # To match the expected directory strcutre, e.g. 15060/ contains depth 1 maps # starting at 15060, make the intermediate directory if not file_path.parent.exists(): file_path.parent.mkdir(parents=True, exist_ok=True) - with open(file_path, "wb") as f: - f.write(response.content) + partial_path = file_path.with_suffix(file_path.suffix + ".part") + try: + with requests.get(url, stream=True, timeout=(30, 120)) as response: + response.raise_for_status() + with partial_path.open("wb") as f: + f.writelines(response.iter_content(1024 * 1024)) + partial_path.replace(file_path) + finally: + partial_path.unlink(missing_ok=True) # Set the location of the file in the cache request.config.cache.set( "downloaded_file_" + os.path.basename(url), str(file_path) ) - return cache_dir + return cache_dir.resolve() def test_act(database_sessionmaker, downloaded_data_file): @@ -240,7 +248,9 @@ def test_sky_coverage(database_sessionmaker, downloaded_data_file): ) act.core(session=database_sessionmaker, args=args) - update_sky_coverage.core(session=database_sessionmaker, convention="ACT") + update_sky_coverage.core( + session=database_sessionmaker, relative_to=args.relative_to + ) with database_sessionmaker() as session: d1maps = session.query(DepthOneMapTable).all() for d1map in d1maps: @@ -276,22 +286,22 @@ def test_sky_coverage_2(database_sessionmaker, downloaded_data_file): ) act.core(session=database_sessionmaker, args=args) - update_sky_coverage.core(session=database_sessionmaker, convention="ACT") + update_sky_coverage.core( + session=database_sessionmaker, relative_to=args.relative_to + ) 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.linspace(0, nonzero_radec.shape[1] - 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 + coord = ICRS(pix[1] * u.rad, pix[0] * u.rad) return_d1map = get_maps_by_coverage(coord, session) assert d1map.map_name in [m.map_name for m in return_d1map] @@ -301,7 +311,7 @@ def test_sky_coverage_2(database_sessionmaker, downloaded_data_file): assert len(return_d1map) == 0 coord2 = ICRS( - 180 * u.rad, 0 * u.rad + 180 * u.deg, 0 * u.deg ) # Test a point on the opposite side of the sky return_d1map_list = get_maps_by_coverage([coord, coord2], session) assert len(return_d1map_list) == 2 @@ -315,7 +325,7 @@ def test_sky_coverage_2(database_sessionmaker, downloaded_data_file): 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 +def test_ra_to_index(): + assert update_sky_coverage.ra_to_index(-180) == 18 + assert update_sky_coverage.ra_to_index(170) == 17 + assert update_sky_coverage.ra_to_index(0) == 0