diff --git a/nuztf/base_scanner.py b/nuztf/base_scanner.py index 6741234..09d0461 100644 --- a/nuztf/base_scanner.py +++ b/nuztf/base_scanner.py @@ -1,6 +1,7 @@ #!/usr/bin/env python3 # coding: utf-8 +import json import logging import time from pathlib import Path @@ -32,7 +33,6 @@ from nuztf.flatpix import get_flatpix, get_nested_pix from nuztf.fritz import save_source_to_group from nuztf.observations import get_obs_summary -from nuztf.observations_depot import get_obs_summary as alt_get_obs_summary from nuztf.paths import BASE_CANDIDATE_DIR, RESULTS_DIR from nuztf.plot import lightcurve_from_alert from nuztf.utils import cosmo @@ -63,7 +63,7 @@ def __init__( self.nside, self.map_probs, self.data, - self.pixel_area, + self.total_pixel_area, self.key, ) = self.unpack_skymap() @@ -120,6 +120,8 @@ def __init__( self.observations = None + self.pixel_area = hp.pixelfunc.nside2pixarea(self.nside, degrees=True) + if not hasattr(self, "dist"): self.dist = None @@ -135,7 +137,7 @@ def get_output_dir(self) -> Path: return output_dir def get_cache_dir(self) -> Path: - cache_dir = BASE_CANDIDATE_DIR.joinpath(self.get_name()) + cache_dir = BASE_CANDIDATE_DIR.joinpath(self.get_name().replace(" ", "_")) cache_dir.mkdir(exist_ok=True, parents=True) return cache_dir @@ -240,9 +242,9 @@ def check_ampel_filter(self, ztf_name): return pipeline_bool def get_multi_night_summary(self, max_days=None): - mns = alt_get_obs_summary(self.t_min, max_days=max_days) - if mns is None: - mns = get_obs_summary(self.t_min, max_days=max_days) + mns = get_obs_summary(self.t_min, max_days=max_days) + # if mns is None: + # mns = get_obs_summary(self.t_min, max_days=max_days) return mns def query_ampel( @@ -327,6 +329,11 @@ def scan_area( """ query_res = self.query_ampel(t_min=t_min, t_max=t_max) + cache_file_initial_stage = self.get_cache_dir().joinpath("initial_stage.json") + + with open(cache_file_initial_stage, "w") as outfile: + json.dump(query_res, outfile) + ztf_ids_first_stage = [] for res in tqdm(query_res): if self.filter_f_no_prv(res): @@ -538,7 +545,7 @@ def extract_npix(nside: int, ra_deg: float, dec_deg: float) -> int: return int(hp.ang2pix(nside, theta, phi, nest=True)) - def create_candidate_summary(self): + def create_candidate_summary(self, include_ps1: bool = True) -> str: """ Create pdf with lightcurve plots of all candidates @@ -560,6 +567,7 @@ def create_candidate_summary(self): fig, _ = lightcurve_from_alert( [alert], include_cutouts=True, + include_ps1=include_ps1, logger=self.logger, t_0_mjd=self.t_min.mjd, ) @@ -791,245 +799,12 @@ def text_summary(self): return text def calculate_overlap_with_observations( - self, fields=None, pid=None, first_det_window_days=3.0, min_sep=0.01 - ): - """ - Calculate the overlap of the candidates with the observations - - :param fields: list of fields to consider. By default, use actual log. - :param pid: pid to consider. By default, use all. - :param first_det_window_days: number of days to consider for first detection - :param min_sep: minimum separation in days between observations - """ - - if fields is None: - mns = self.get_multi_night_summary(first_det_window_days) - - else: - - class MNS: - def __init__(self, data): - self.data = pandas.DataFrame( - data, columns=["field", "ra", "dec", "datetime"] - ) - - data = [] - - for f in fields: - ra, dec = ztfquery_fields.get_field_centroid(f)[0] - for i in range(2): - t = Time(self.t_min.jd + 0.1 * i, format="jd").utc - t.format = "isot" - t = t.value - data.append([f, ra, dec, t]) - - mns = MNS(data) - - # Skip all 64 simultaneous quadrant entries, we only need one per observation for qa log - - data = mns.data.copy() - - ras = np.ones_like(data["field"]) * np.nan - decs = np.ones_like(data["field"]) * np.nan - - # Actually load up RA/Dec - - veto_fields = [] - - for field in list(set(data["field"])): - mask = data["field"] == field - - res = ztfquery_fields.get_field_centroid(field) - - if len(res) > 0: - ras[mask] = res[0][0] - decs[mask] = res[0][1] - - else: - veto_fields.append(field) - - if len(veto_fields) > 0: - self.logger.info( - f"No RA/Dec found by ztfquery for fields {veto_fields}. " - f"These observation have to be ignored." - ) - - data["ra"] = ras - data["dec"] = decs - - mask = np.array([~np.isnan(x) for x in data["ra"]]) - - data = data[mask] - - if pid is not None: - pid_mask = data["pid"] == str(pid) - data = data[pid_mask] - - obs_times = np.array( - [ - Time( - data["datetime"].iat[i].replace(" ", "T"), - format="isot", - scale="utc", - ) - for i in range(len(data)) - ] - ) - - if first_det_window_days is not None: - first_det_mask = [ - x < Time(self.t_min.jd + first_det_window_days, format="jd").utc - for x in obs_times - ] - data = data[first_det_mask] - obs_times = obs_times[first_det_mask] - - if len(obs_times) == 0: - self.logger.warning("No observations found for this event.") - return None - - self.logger.info(f"Most recent observation found is {obs_times[-1]}") - self.logger.info("Unpacking observations") - - # re-read the map and subsample it to nside=64 if nside is too large, as overlap calculation for big maps with large nside take forever - if self.nside > 256: - ( - self.map_coords, - self.pixel_nos, - self.nside, - self.map_probs, - self.data, - self.pixel_area, - self.key, - ) = self.unpack_skymap(output_nside=256) - - pix_map = dict() - pix_obs_times = dict() - - field_pix = get_flatpix(nside=self.nside, logger=self.logger) - - for i, obs_time in enumerate(tqdm(obs_times)): - field = data["field"].iat[i] - - flat_pix = field_pix[field] - - t = obs_time.jd - - for p in flat_pix: - if p not in pix_obs_times.keys(): - pix_obs_times[p] = [t] - else: - pix_obs_times[p] += [t] - - if p not in pix_map.keys(): - pix_map[p] = [field] - else: - pix_map[p] += [field] - - npix = hp.nside2npix(self.nside) - theta, phi = hp.pix2ang(self.nside, np.arange(npix), nest=False) - radecs = SkyCoord(ra=phi * u.rad, dec=(0.5 * np.pi - theta) * u.rad) - idx = np.where(np.abs(radecs.galactic.b.deg) <= 10.0)[0] - - double_in_plane_pixels = [] - double_in_plane_probs = [] - single_in_plane_pixels = [] - single_in_plane_prob = [] - veto_pixels = [] - plane_pixels = [] - plane_probs = [] - times = [] - double_no_plane_prob = [] - double_no_plane_pixels = [] - single_no_plane_prob = [] - single_no_plane_pixels = [] - - overlapping_fields = [] - - for i, p in enumerate(tqdm(hp.nest2ring(self.nside, self.pixel_nos))): - if p in pix_obs_times.keys(): - if p in idx: - plane_pixels.append(p) - plane_probs.append(self.map_probs[i]) - - obs = pix_obs_times[p] - - # check which healpix are observed twice - if max(obs) - min(obs) > min_sep: - # is it in galactic plane or not? - if p not in idx: - double_no_plane_prob.append(self.map_probs[i]) - double_no_plane_pixels.append(p) - else: - double_in_plane_probs.append(self.map_probs[i]) - double_in_plane_pixels.append(p) - - else: - if p not in idx: - single_no_plane_pixels.append(p) - single_no_plane_prob.append(self.map_probs[i]) - else: - single_in_plane_prob.append(self.map_probs[i]) - single_in_plane_pixels.append(p) - - overlapping_fields += pix_map[p] - - times += list(obs) - else: - veto_pixels.append(p) - - overlapping_fields = sorted(list(set(overlapping_fields))) - - _observations = data.query("obsjd in @times").reset_index(drop=True)[ - ["obsjd", "datetime", "exptime", "fid", "field"] - ] - bands = [self.fid_to_band(fid) for fid in _observations["fid"].values] - _observations["band"] = bands - _observations.drop(columns=["fid"], inplace=True) - self.observations = _observations - - self.logger.info("All observations:") - self.logger.info(f"\n{self.observations}") - - try: - self.first_obs = Time(min(times), format="jd") - self.first_obs.utc.format = "isot" - self.last_obs = Time(max(times), format="jd") - self.last_obs.utc.format = "isot" - - except ValueError: - err = ( - f"No observations of this field were found at any time between {self.t_min} and" - f"{obs_times[-1]}. Coverage overlap is 0%, but recent observations might be missing!" - ) - self.logger.error(err) - raise ValueError(err) - - self.logger.info(f"Observations started at {self.first_obs.jd}") - - return ( - double_in_plane_pixels, - double_in_plane_probs, - single_in_plane_pixels, - single_in_plane_prob, - veto_pixels, - plane_pixels, - plane_probs, - times, - double_no_plane_prob, - double_no_plane_pixels, - single_no_plane_prob, - single_no_plane_pixels, - overlapping_fields, - ) - - def calculate_overlap_with_depot_observations( self, first_det_window_days=3.0, min_sep=0.01 ): - mns = alt_get_obs_summary(t_min=self.t_min, max_days=first_det_window_days) + mns = get_obs_summary(t_min=self.t_min, max_days=first_det_window_days) if mns is None: - return None + return None, None, None data = mns.data.copy() @@ -1048,7 +823,7 @@ def calculate_overlap_with_depot_observations( self.nside, self.map_probs, self.data, - self.pixel_area, + self.total_pixel_area, self.key, ) = self.unpack_skymap(output_nside=256) @@ -1063,13 +838,13 @@ def calculate_overlap_with_depot_observations( field = obs["field_id"].iloc[0] try: - flat_pix = nested_pix[field] + flat_pix = nested_pix[int(field)] - mask = obs["status"] == 0 + mask = obs["status"].astype(int) == 0 indices = obs["qid"].values[mask] for qid in indices: - pixels = flat_pix[qid] + pixels = flat_pix[int(qid)] for p in pixels: if p not in pix_obs_times.keys(): @@ -1082,7 +857,7 @@ def calculate_overlap_with_depot_observations( else: pix_map[p] += [field] - except KeyError: + except (KeyError, ValueError): self.logger.warning( f"Field {field} not found in nested pix dict. " f"This might be an engineering observation." @@ -1091,54 +866,43 @@ def calculate_overlap_with_depot_observations( npix = hp.nside2npix(self.nside) theta, phi = hp.pix2ang(self.nside, np.arange(npix), nest=False) radecs = SkyCoord(ra=phi * u.rad, dec=(0.5 * np.pi - theta) * u.rad) - idx = np.where(np.abs(radecs.galactic.b.deg) <= 10.0)[0] - - double_in_plane_pixels = [] - double_in_plane_probs = [] - single_in_plane_pixels = [] - single_in_plane_prob = [] - veto_pixels = [] - plane_pixels = [] - plane_probs = [] - times = [] - double_no_plane_prob = [] - double_no_plane_pixels = [] - single_no_plane_prob = [] - single_no_plane_pixels = [] + + ras, decs = hp.pixelfunc.pix2ang( + self.nside, hp.nest2ring(self.nside, self.pixel_nos), lonlat=True + ) + + coverage_data = [] overlapping_fields = [] + times = [] for i, p in enumerate(tqdm(hp.nest2ring(self.nside, self.pixel_nos))): - if p in pix_obs_times.keys(): - if p in idx: - plane_pixels.append(p) - plane_probs.append(self.map_probs[i]) + entry = { + "pixel_no": p, + "prob": self.map_probs[i], + "gal_b": radecs[i].galactic.b.deg, + "in_plane": abs(radecs[i].galactic.b.deg) < 10.0, + "ra_deg": ras[i], + "dec_deg": decs[i], + } + if p in pix_obs_times.keys(): obs = pix_obs_times[p] - # check which healpix are observed twice - if max(obs) - min(obs) > min_sep: - # is it in galactic plane or not? - if p not in idx: - double_no_plane_prob.append(self.map_probs[i]) - double_no_plane_pixels.append(p) - else: - double_in_plane_probs.append(self.map_probs[i]) - double_in_plane_pixels.append(p) - - else: - if p not in idx: - single_no_plane_pixels.append(p) - single_no_plane_prob.append(self.map_probs[i]) - else: - single_in_plane_prob.append(self.map_probs[i]) - single_in_plane_pixels.append(p) + entry["n_det_class"] = 2 if max(obs) - min(obs) > min_sep else 1 + + entry["latency_days"] = min(obs) - self.t_min.jd overlapping_fields += pix_map[p] times += list(obs) else: - veto_pixels.append(p) + entry["n_det_class"] = 0 + entry["latency_days"] = np.nan + + coverage_data.append(entry) + + coverage_df = pd.DataFrame(coverage_data) overlapping_fields = sorted(list(set(overlapping_fields))) @@ -1171,18 +935,8 @@ def calculate_overlap_with_depot_observations( self.logger.info(f"Observations started at {self.first_obs.isot}") return ( - double_in_plane_pixels, - double_in_plane_probs, - single_in_plane_pixels, - single_in_plane_prob, - veto_pixels, - plane_pixels, - plane_probs, + coverage_df, times, - double_no_plane_prob, - double_no_plane_pixels, - single_no_plane_prob, - single_no_plane_pixels, overlapping_fields, ) @@ -1195,100 +949,75 @@ def plot_overlap_with_observations(self, first_det_window_days=None, min_sep=0.0 """ - overlap_res = self.calculate_overlap_with_depot_observations( + ( + coverage_df, + times, + overlapping_fields, + ) = self.calculate_overlap_with_observations( first_det_window_days=first_det_window_days, min_sep=min_sep, ) - if overlap_res is None: - self.logger.info("IPAC depot failed, using ztfquery to obtain observations") - overlap_res = self.calculate_overlap_with_observations( - first_det_window_days=first_det_window_days, - min_sep=min_sep, - ) - if overlap_res is None: + if coverage_df is None: self.logger.warning("Not plotting overlap with observations.") return - else: - ( - double_in_plane_pixels, - double_in_plane_probs, - single_in_plane_pixels, - single_in_plane_prob, - veto_pixels, - plane_pixels, - plane_probs, - times, - double_no_plane_prob, - double_no_plane_pixels, - single_no_plane_prob, - single_no_plane_pixels, - overlapping_fields, - ) = overlap_res fig = plt.figure() plt.subplot(projection="aitoff") self.overlap_fields = list(set(overlapping_fields)) - self.overlap_prob = np.sum(double_in_plane_probs + double_no_plane_prob) * 100.0 + overlap_mask = (coverage_df["n_det_class"] == 2) & ~coverage_df["in_plane"] + self.overlap_prob = coverage_df[overlap_mask]["prob"].sum() * 100.0 size = hp.max_pixrad(self.nside) ** 2 * 50.0 - veto_pos = np.array( - [hp.pixelfunc.pix2ang(self.nside, i, lonlat=True) for i in veto_pixels] - ).T + veto_pixels = coverage_df.where(coverage_df["n_det_class"] == 0) - if len(veto_pos) > 0: + if len(veto_pixels) > 0: plt.scatter( - np.radians(self.wrap_around_180(veto_pos[0])), - np.radians(veto_pos[1]), + np.radians(veto_pixels["ra_deg"]), + np.radians(veto_pixels["dec_deg"]), color="red", s=size, ) - plane_pos = np.array( - [hp.pixelfunc.pix2ang(self.nside, i, lonlat=True) for i in plane_pixels] - ).T + plane_pixels = coverage_df.where( + coverage_df["in_plane"] & (coverage_df["n_det_class"] > 0) + ) - if len(plane_pos) > 0: + if len(plane_pixels) > 0: plt.scatter( - np.radians(self.wrap_around_180(plane_pos[0])), - np.radians(plane_pos[1]), + np.radians(self.wrap_around_180(plane_pixels["ra_deg"])), + np.radians(plane_pixels["dec_deg"]), color="green", s=size, ) - single_pos = np.array( - [ - hp.pixelfunc.pix2ang(self.nside, i, lonlat=True) - for i in single_no_plane_pixels - ] - ).T + single_pixels = coverage_df.where( + ~coverage_df["in_plane"] & (coverage_df["n_det_class"] == 1) + ) - if len(single_pos) > 0: + if len(single_pixels) > 0: plt.scatter( - np.radians(self.wrap_around_180(single_pos[0])), - np.radians(single_pos[1]), - c=single_no_plane_prob, + np.radians(self.wrap_around_180(single_pixels["ra_deg"])), + np.radians(single_pixels["dec_deg"]), + c=single_pixels["prob"], vmin=0.0, vmax=max(self.data[self.key]), s=size, cmap="gray", ) - plot_pos = np.array( - [ - hp.pixelfunc.pix2ang(self.nside, i, lonlat=True) - for i in double_no_plane_pixels - ] - ).T + double_pixels = coverage_df.where( + ~coverage_df["in_plane"] & (coverage_df["n_det_class"] == 2) + ) - if len(plot_pos) > 0: + if len(double_pixels) > 0: plt.scatter( - np.radians(self.wrap_around_180(plot_pos[0])), - np.radians(plot_pos[1]), - c=double_no_plane_prob, + np.radians(self.wrap_around_180(double_pixels["ra_deg"])), + np.radians(double_pixels["dec_deg"]), + c=double_pixels["prob"], vmin=0.0, vmax=max(self.data[self.key]), s=size, @@ -1309,35 +1038,20 @@ def plot_overlap_with_observations(self, first_det_window_days=None, min_sep=0.0 "In total, {3:.1f} % of the contour was observed at least twice, " "and excluding low galactic latitudes.\n" "These estimates account for chip gaps.".format( - 100 - * ( - np.sum(double_in_plane_probs) - + np.sum(single_in_plane_prob) - + np.sum(single_no_plane_prob) - + np.sum(double_no_plane_prob) - ), - 100 * np.sum(plane_probs), - 100.0 * (np.sum(double_in_plane_probs) + np.sum(double_no_plane_prob)), - 100.0 * np.sum(double_no_plane_prob), + 100.0 * coverage_df.query("n_det_class > 0")["prob"].sum(), + 100.0 * coverage_df.query("in_plane & n_det_class > 0")["prob"].sum(), + 100.0 * coverage_df.query("n_det_class == 2")["prob"].sum(), + 100.0 * coverage_df.query("n_det_class == 2 & ~in_plane")["prob"].sum(), ) ) - n_pixels = len( - single_in_plane_pixels - + double_in_plane_pixels - + double_no_plane_pixels - + single_no_plane_pixels - ) - n_double = len(double_no_plane_pixels + double_in_plane_pixels) - n_plane = len(plane_pixels) + n_pixels = len(coverage_df.query("n_det_class > 0 & prob > 0.0")) + n_double = len(coverage_df.query("n_det_class == 2 & prob > 0.0")) + n_plane = len(coverage_df.query("in_plane & n_det_class > 0 & prob > 0.0")) - self.healpix_area = ( - hp.pixelfunc.nside2pixarea(self.nside, degrees=True) * n_pixels - ) - self.double_extragalactic_area = ( - hp.pixelfunc.nside2pixarea(self.nside, degrees=True) * n_double - ) - plane_area = hp.pixelfunc.nside2pixarea(self.nside, degrees=True) * n_plane + self.healpix_area = self.pixel_area * n_pixels + self.double_extragalactic_area = self.pixel_area * n_double + plane_area = self.pixel_area * n_plane self.overlap_fields = overlapping_fields @@ -1367,7 +1081,7 @@ def get_exposure_summary(self) -> pd.DataFrame: max_days = max(self.observations["obsjd"]) - self.t_min.jd - mns = alt_get_obs_summary(t_min=self.t_min, max_days=max_days) + mns = get_obs_summary(t_min=self.t_min, max_days=max_days) exposures = [] diff --git a/nuztf/irsa.py b/nuztf/irsa.py index 2e09e14..b7bd2d6 100644 --- a/nuztf/irsa.py +++ b/nuztf/irsa.py @@ -18,7 +18,6 @@ from ztfquery.lightcurve import LCQuery from nuztf.ampel_api import ampel_api_name -from nuztf.observations import get_most_recent_obs from nuztf.parse_nu_gcn import find_gcn_no, parse_gcn_circular from nuztf.style import base_height, base_width, big_fontsize, dpi, plot_dir from nuztf.utils import cosmo, is_tns_name, is_ztf_name, query_tns_by_name @@ -307,24 +306,6 @@ def plot_irsa_lightcurve( f"{latest['filtercode'][1]}={latest['mag']:.2f}+/-{latest['magerr']:.2f}" ) - # If you want, you can check the most recent observation - - if check_obs: - logger.info( - f"Getting most recent ZTF observation, looking back {check_obs_lookback_weeks} weeks." - ) - mro = get_most_recent_obs( - ra=source_coords[0], - dec=source_coords[1], - lookback_weeks_max=check_obs_lookback_weeks, - ) - - if mro is not None: - ot = format_date(Time(mro["obsjd"], format="jd"), atel=atel) - logger.info(f"Most recent observation at {ot}") - else: - logger.info("No recent observation found.") - # Plot each band (g/r/i) for fc in ["zg", "zr", "zi"]: diff --git a/nuztf/neutrino_scanner.py b/nuztf/neutrino_scanner.py index 14ae8db..425d2ce 100644 --- a/nuztf/neutrino_scanner.py +++ b/nuztf/neutrino_scanner.py @@ -250,8 +250,16 @@ def unpack_skymap(self, skymap=None, output_nside: None | int = None): data = np.zeros(hp.nside2npix(output_nside), dtype=np.dtype([(key, float)])) data[np.array(pixel_nos)] = map_probs - pixel_area = hp.nside2pixarea(output_nside, degrees=True) * float( + total_pixel_area = hp.nside2pixarea(output_nside, degrees=True) * float( len(map_coords) ) - return map_coords, pixel_nos, output_nside, map_probs, data, pixel_area, key + return ( + map_coords, + pixel_nos, + output_nside, + map_probs, + data, + total_pixel_area, + key, + ) diff --git a/nuztf/observations.py b/nuztf/observations.py index 54a1ace..0084c3d 100644 --- a/nuztf/observations.py +++ b/nuztf/observations.py @@ -1,10 +1,9 @@ -import glob +import json import logging import os import time -import warnings from glob import glob -from typing import Optional +from pathlib import Path import backoff import numpy as np @@ -14,188 +13,381 @@ from astropy import units as u from astropy.time import Time from pyvo.auth import authsession, securitymethods +from requests.auth import HTTPBasicAuth +from requests.exceptions import HTTPError +from tqdm import tqdm from ztfquery import skyvision -from ztfquery.fields import get_fields_containing_target from ztfquery.io import LOCALSOURCE from nuztf import credentials -from nuztf.fritz import fritz_api -logger = logging.getLogger(__name__) +coverage_dir = Path(LOCALSOURCE).joinpath("all_obs") +coverage_dir.mkdir(exist_ok=True) -username, password = credentials.load_credentials("irsa") +partial_flag = "_PARTIAL" -# Gather login information -data = {"username": username, "password": password} +logger = logging.getLogger(__name__) +# IPAC TAP login +username, password = credentials.load_credentials("ipacdepot") +data = {"username": username, "password": password} headers = {"Content-Type": "application/x-www-form-urlencoded", "Accept": "text/plain"} +IPAC_TAP_URL = "https://irsa.ipac.caltech.edu" -login_url = "https://irsa.ipac.caltech.edu" -coverage_dir = os.path.join(LOCALSOURCE, "all_obs") +class MNS: + def __init__(self, df): + self.data = df -try: - os.makedirs(coverage_dir) -except OSError: - pass -partial_flag = "_PARTIAL" +class NoDepotEntry(Exception): + """ + No entry in the depot for a given date + """ + + +def get_date(jd: float) -> str: + """ + Get the date in YYYYMMDD format from a JD + + :param jd: JD + :return: String date + """ + return str(Time(jd, format="jd").isot).split("T")[0].replace("-", "") + + +def download_depot_log(date): + """ + Download the depot log for a given date + + :param date: date in YYYYMMDD format + :return: json log + """ + url = f"https://ztfweb.ipac.caltech.edu/ztf/depot/{date}/ztf_recentproc_{date}.json" + response = requests.get(url, auth=HTTPBasicAuth(username, password)) + if response.status_code == 404: + raise NoDepotEntry(f"No depot entry for {date} at url {url}") + response.raise_for_status() + return response.json() -def coverage_path(jd: float) -> str: +# The following functions are used to find the paths to the cached coverage files + + +def coverage_depot_path(jd: float) -> Path: + """ + Return the path to the cached coverage file for a given JD + + :param jd: JD + :return: path to cached coverage file + """ + + date = get_date(jd) + if (Time.now().jd - jd) < 1: partial_ext = partial_flag else: partial_ext = "" - return os.path.join(coverage_dir, f"{jd}{partial_ext}.csv") + return coverage_dir.joinpath(f"{date}{partial_ext}.json") -class MNS: - def __init__(self, df): - self.data = df +def coverage_skyvision_path(jd: float) -> Path: + """ + Find the path to the Skyvision coverage file for a given JD + + :param jd: JD + :return: Output path + """ + if (Time.now().jd - jd) < 1: + partial_ext = partial_flag + else: + partial_ext = "" + + output_path = coverage_dir.joinpath(f"{get_date(jd)}{partial_ext}_skyvision.json") + + return output_path + + +def coverage_tap_path(jd: float) -> Path: + """ + Find the path to the TAP coverage file for a given JD + + :param jd: JD + :return: Output path + """ + if (Time.now().jd - jd) < 1: + partial_ext = partial_flag + else: + partial_ext = "" + + output_path = coverage_dir.joinpath(f"{get_date(jd)}{partial_ext}_tap.json") + + return output_path + + +# The following functions are used to write the coverage to the cache + + +def write_coverage_depot(jds: list[float]): + """ + Write the depot coverage for a list of JDs to the cache + + Requires a valid IPAC depot login + + :param jds: JDs + :return: None + """ + for jd in tqdm(jds): + date = get_date(jd) + try: + log = download_depot_log(date) + except NoDepotEntry as exc: + log = {} + logger.warning(f"No depot entry for {date}: {exc}") + + path = coverage_depot_path(jd) + with open(path, "w") as f: + json.dump(log, f) @backoff.on_exception(backoff.expo, requests.exceptions.RequestException, max_time=60) -def write_coverage(jds: [int]): +def write_coverage_tap(jds: [float]): + """ + Write the coverage for a list of JDs to the cache using TAP + + (JDs must be half-integer) + + :param jds: Half-integer JDs + :return: None + """ with requests.Session() as session: # Create a session and do the login. # The cookie will end up in the cookie jar of the session. - response = session.post(login_url, data=data, headers=headers) + response = session.post(IPAC_TAP_URL, data=data, headers=headers) response.raise_for_status() auth = authsession.AuthSession() auth.credentials.set(securitymethods.ANONYMOUS, session) client = pyvo.dal.TAPService("https://irsa.ipac.caltech.edu/TAP", auth) for jd in jds: + assert jd - int(jd) == 0.5, "JD must be a half-integer" + obstable = client.search( f""" - SELECT field,rcid,fid,expid,obsjd,exptime,maglimit,ipac_gid,seeing + SELECT expid,obsjd,fid,field,exptime,rcid,maglimit,infobits,ipac_gid FROM ztf.ztf_current_meta_sci WHERE (obsjd BETWEEN {jd} AND {jd + 1}) """ ).to_table() - names = ("ipac_gid",) - renames = ("programid",) + names = ( + "expid", + "exptime", + "field", + "fid", + "rcid", + "maglimit", + "infobits", + "ipac_gid", + ) + renames = ( + "exposure_id", + "exposure_time", + "field_id", + "filter_id", + "qid", + "maglim", + "status", + "programid", + ) obstable.rename_columns(names, renames) - obs_grouped_by_exp = obstable.to_pandas().groupby("expid") + obs = obstable.to_pandas() + obs["nalertpackets"] = np.nan - output_path = coverage_path(jd) + output_path = coverage_tap_path(jd) - with open(output_path, "w") as fid: - fid.write( - "obsid,field,obsjd,datetime,seeing,limmag,exposure_time,fid,processed_fraction\n" - ) + obs.to_json(output_path) + time.sleep(1) - for group_name, df_group in obs_grouped_by_exp: - processed_fraction = len(df_group["field"]) / 64.0 + logger.warning( + f"Coverage for JD {jd} written to {output_path} using TAP, " + f"but this will only include programid=1" + ) - t = Time(df_group["obsjd"].iloc[0], format="jd") - t.format = "isot" - line = f'{int(df_group["expid"].iloc[0])},{int(df_group["field"].iloc[0])},{t.jd},{t},{df_group["seeing"].median():.10f},{df_group["maglimit"].median():.10f},{df_group["exptime"].iloc[0]},{int(df_group["fid"].iloc[0])},{processed_fraction} \n' - fid.write(line) +def write_coverage_skyvision(jds: list[float]): + """ + Write the Skyvision coverage for a list of JDs to the cache - time.sleep(1) + :param jds: JDs + :return: None + """ + for jd in tqdm(jds): + date = get_date(jd) + + assert jd - int(jd) == 0.5, "JD must be a half-integer" + + res = skyvision.get_log(f"{date[:4]}-{date[4:6]}-{date[6:8]}", verbose=False) + + path = coverage_skyvision_path(jd) + + if res is not None: + # The skyvision log has a bug where some entries have a FieldID of "NONE" + mask = ( + pd.notnull(res["FieldID"]) & res["FieldID"].astype(str).str.isnumeric() + ) + res = res[mask] + jds = [ + Time(f'{row["UT Date"]}T{row["UT Time"]}').jd + for _, row in res.iterrows() + ] + res["obsjd"] = jds + res["status"] = (res["Observation Status"] == "FAILED").astype(int) + res["filter_id"] = res["Filter"].apply( + lambda x: 1 if x == "FILTER_ZTF_G" else 2 if x == "FILTER_ZTF_R" else 3 + ) + res["maglim"] = 20.5 + res["field_id"] = res["FieldID"].astype(int) + res["exposure_time"] = res["Exptime"] + res = res[ + ["obsjd", "filter_id", "field_id", "exposure_time", "maglim", "status"] + ] + + new_res = [] + for _, row in res.iterrows(): + new = row.to_dict() + new_res += [dict(qid=int(i), **new) for i in range(64)] + pd.DataFrame(new_res).to_json(path) + + else: + with open(path, "w") as f: + json.dump({}, f) + + +# The following functions are used to get the coverage from the cache -def get_coverage(jds: [int]) -> Optional[pd.DataFrame]: +def get_coverage(jds: [int]) -> pd.DataFrame | None: + """ + Get a dataframe of the coverage for a list of JDs + + Will use the cache if available, otherwise will query the depot, and lastly TAP + + :param jds: JDs + :return: Coverage dataframe + """ # Clear any logs flagged as partial/incomplete - cache_files = glob(f"{coverage_dir}/*.csv") - partial_logs = [x for x in cache_files if partial_flag in x] + cache_files = coverage_dir.glob("*.json") + partial_logs = [x for x in cache_files if partial_flag in str(x)] if len(partial_logs) > 0: logger.debug(f"Removing the following partial logs: {partial_logs}") for partial_log in partial_logs: - os.remove(partial_log) + partial_log.unlink() # Only write missing logs missing_logs = [] for jd in jds: - if not os.path.exists(coverage_path(jd)): + depot_path = coverage_depot_path(jd) + if not depot_path.exists(): missing_logs.append(jd) else: - df = pd.read_csv(coverage_path(jd)) + df = pd.read_json(coverage_depot_path(jd)) if len(df) == 0: missing_logs.append(jd) if len(missing_logs) > 0: - logger.debug( - f"Some logs were missing from the cache. Querying for the following JDs: {missing_logs}" + logger.info( + f"Some logs were missing from the cache. " + f"Querying for the following JDs in depot: {missing_logs}" ) - write_coverage(missing_logs) + write_coverage_depot(missing_logs) + + # Try skyvision for missing logs + still_missing_logs = [] + for jd in missing_logs: + skyvision_path = coverage_skyvision_path(jd) + if not skyvision_path.exists(): + still_missing_logs.append(jd) + else: + df = pd.read_json(coverage_skyvision_path(jd)) + if len(df) == 0: + still_missing_logs.append(jd) + + if len(still_missing_logs) > 0: + logger.info( + f"Some logs were still missing from the cache. " + f"Querying for the following JDs from skyvision: {still_missing_logs}" + ) + write_coverage_skyvision(still_missing_logs) + + # Try TAP for missing logs + + completely_missing_logs = [] + + for jd in still_missing_logs: + depot_path = coverage_depot_path(jd) + if not depot_path.exists(): + completely_missing_logs.append(jd) + else: + df = pd.read_json(coverage_depot_path(jd)) + if len(df) == 0: + completely_missing_logs.append(jd) + + if len(completely_missing_logs) > 0: + logger.info( + f"Some logs were still missing from the cache. " + f"Querying for the following JDs from TAP: {completely_missing_logs}" + ) + write_coverage_tap(completely_missing_logs) # Load logs from cache results = [] - for jd in jds: - res = pd.read_csv(coverage_path(jd)) - results.append(res) + for jd in tqdm(jds): + res = pd.read_json(coverage_depot_path(jd)) + if len(res) > 0: + results.append(res) + else: + res = pd.read_json(coverage_skyvision_path(jd)) + if len(res) > 0: + results.append(res) + else: + res = pd.read_json(coverage_tap_path(jd)) + results.append(res) if results: - result_df = pd.concat(results) + result_df = pd.concat(results, ignore_index=True) return result_df else: return None -def get_obs_summary(t_min, t_max=None, max_days: int = None): +def get_obs_summary_depot(t_min: Time, t_max: Time) -> MNS | None: """ - Get observation summary from Skyvision (or IRSA if Skyvision fails) + Get observation summary from depot """ - now = Time.now() - if t_max and max_days: - raise ValueError("Choose either t_max or max_days, not both") + jds = np.arange(int(t_min.jd) - 0.5, int(t_max.jd) + 1.5) - if t_max is None: - if max_days is None: - t_max = now - else: - t_max = t_min + (max_days * u.day) - - if t_max > now: - t_max = now - - logger.info("Getting observation logs from skyvision.") - mns = get_obs_summary_skyvision(t_min=t_min, t_max=t_max) - - # Useless for now as long as Fritz also depends on TAP - - # if len(mns.data) == 0: - # logger.debug("Empty observation log, try Fritz instead.") - # mns = get_obs_summary_fritz(t_min=t_min, t_max=t_max) + res = get_coverage(jds) - if len(mns.data) == 0: - logger.debug("Empty observation log, try IRSA instead.") - mns = get_obs_summary_irsa(t_min=t_min, t_max=t_max) - - logger.debug(f"Found {len(mns.data)} observations in total.") - - return mns - - -def get_obs_summary_irsa(t_min, t_max): - """ - Get observation summary from IRSA - """ - jds = np.arange(int(t_min.jd), int(t_max.jd) + 1) + if len(res) == 0: + return None - logger.debug("Getting coverage") + res["date"] = Time(res["obsjd"].to_numpy(), format="jd").isot - df = get_coverage(jds) - mns = MNS(df) + mns = MNS(df=res) mns.data.query(f"obsjd >= {t_min.jd} and obsjd <= {t_max.jd}", inplace=True) + mns.data.reset_index(inplace=True) mns.data.drop(columns=["index"], inplace=True) - logger.debug("Done") - return mns @@ -237,101 +429,30 @@ def get_obs_summary_skyvision(t_min, t_max): return mns -def get_obs_summary_fritz(t_min, t_max) -> MNS | None: +def get_obs_summary(t_min, t_max=None, max_days: float = None) -> MNS | None: """ - Get observation summary from Fritz + Get observation summary from IPAC depot """ + now = Time.now() - res = fritz_api(method="GET", endpoint_extension="api/observation", data=data) - res = res.json().get("data", {}).get("observations") - - t_min_date = t_min.to_value("iso").split(".")[0] - t_max_date = t_max.to_value("iso").split(".")[0] - - logger.debug( - f"Fritz: Obtaining nightly observation logs from {t_min_date} to {t_max_date}" - ) - - params = { - "telescopeName": "Palomar 1.2m Oschin", - "startDate": t_min_date, - "endDate": t_max_date, - } - - res = fritz_api(method="GET", endpoint_extension="api/observation", data=params) - res = res.json().get("data", {}).get("observations") - - if res is None: - return res - - resdict = {} - for i, entry in enumerate(res): - resdict.update( - { - i: { - "datetime": entry["obstime"], - "date": Time(entry["obstime"]).to_value("iso", subfmt="date"), - "exptime": entry["exposure_time"], - "fid": entry["id"], - "field": entry["field"]["field_id"], - "obsjd": Time(entry["obstime"]).to_value("jd"), - "maglim": entry["limmag"], - } - } - ) - - df = pd.DataFrame.from_dict(resdict, orient="index") - - mns = MNS(df) - - logger.debug(f"Found {len(mns.data)} observations in total.") - - return mns - - -def get_most_recent_obs(ra: float, dec: float, lookback_weeks_max: int = 12): - """ """ - - fields = get_fields_containing_target(ra, dec)._data - - logger.info(f"Target in fields {fields}") - - mask = 0.0 - - t_max = Time.now() - - lookback_weeks = 1 - - while np.sum(mask) < 1 and lookback_weeks <= lookback_weeks_max: - t_min = t_max - 1 * u.week + if t_max and max_days: + raise ValueError("Choose either t_max or max_days, not both") - if lookback_weeks > 1: - logger.debug( - f"Searching for most recent obs during the last {lookback_weeks_max} weeks. Looking back {lookback_weeks} weeks." - ) + if t_max is None: + if max_days is None: + t_max = now else: - logger.debug( - f"Searching for most recent obs during the last {lookback_weeks_max} weeks. Looking back {lookback_weeks} week." - ) - - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", category=UserWarning) - - mns = get_obs_summary(t_min=t_min, t_max=t_max) - - mask = np.array([x in fields for x in mns.data["field"]]) - - t_max = t_min + t_max = t_min + (max_days * u.day) - lookback_weeks += 1 + if t_max > now: + t_max = now - if np.sum(mask) >= 1: - index = list(mns.data["datetime"]).index(max(mns.data["datetime"][mask])) - mro = mns.data.iloc[index] - return mro + logger.info("Getting observation logs from IPAC depot.") + mns = get_obs_summary_depot(t_min=t_min, t_max=t_max) + if mns is not None: + logger.info(f"Found {len(set(mns.data))} observations in total.") else: - logger.warning( - f"Found no observation during the last {lookback_weeks_max} weeks." - ) - return None + logger.warning("Found no observations on IPAC depot or TAP.") + + return mns diff --git a/nuztf/observations_depot.py b/nuztf/observations_depot.py deleted file mode 100644 index 2a5c12d..0000000 --- a/nuztf/observations_depot.py +++ /dev/null @@ -1,180 +0,0 @@ -import json -import logging -import os -from glob import glob - -import numpy as np -import pandas as pd -import requests -from astropy import units as u -from astropy.time import Time -from requests.auth import HTTPBasicAuth -from requests.exceptions import HTTPError -from tqdm import tqdm - -from nuztf import credentials -from nuztf.observations import MNS, coverage_dir, partial_flag - -logger = logging.getLogger(__name__) - -username, password = credentials.load_credentials("ipacdepot") - - -class NoDepotEntry(Exception): - """ - No entry in the depot for a given date - """ - - -def download_depot_log(date): - """ - Download the depot log for a given date - - :param date: date in YYYYMMDD format - :return: json log - """ - url = f"https://ztfweb.ipac.caltech.edu/ztf/depot/{date}/ztf_recentproc_{date}.json" - response = requests.get(url, auth=HTTPBasicAuth(username, password)) - if response.status_code == 404: - raise NoDepotEntry(f"No depot entry for {date}") - response.raise_for_status() - return response.json() - - -def coverage_depot_cache(jd: float) -> str: - """ - Return the path to the cached coverage file for a given JD - - :param jd: JD - :return: path to cached coverage file - """ - if (Time.now().jd - jd) < 1: - partial_ext = partial_flag - else: - partial_ext = "" - - return os.path.join(coverage_dir, f"{jd}{partial_ext}.json") - - -def write_depot_coverage(jds: list[int]): - """ - Write the depot coverage for a list of JDs to the cache - - :param jds: JDs - :return: None - """ - for jd in tqdm(jds): - try: - date = str(Time(jd, format="jd").isot).split("T")[0].replace("-", "") - log = download_depot_log(date) - except NoDepotEntry: - log = {} - path = coverage_depot_cache(jd) - with open(path, "w") as f: - json.dump(log, f) - - -def get_coverage_depot(jds: [int]) -> pd.DataFrame | None: - """ - Get a dataframe of the depot coverage for a list of JDs - - :param jds: JDs - :return: Coverage dataframe - """ - # Clear any logs flagged as partial/incomplete - - cache_files = glob(f"{coverage_dir}/*.json") - partial_logs = [x for x in cache_files if partial_flag in x] - - if len(partial_logs) > 0: - logger.debug(f"Removing the following partial logs: {partial_logs}") - for partial_log in partial_logs: - os.remove(partial_log) - - # Only write missing logs - - missing_logs = [] - - for jd in jds: - if not os.path.exists(coverage_depot_cache(jd)): - missing_logs.append(jd) - else: - df = pd.read_json(coverage_depot_cache(jd)) - if len(df) == 0: - missing_logs.append(jd) - - if len(missing_logs) > 0: - logger.debug( - f"Some logs were missing from the cache. " - f"Querying for the following JDs: {missing_logs}" - ) - write_depot_coverage(missing_logs) - - # Load logs from cache - - results = [] - - for jd in tqdm(jds): - res = pd.read_json(coverage_depot_cache(jd)) - results.append(res) - - if results: - result_df = pd.concat(results, ignore_index=True) - return result_df - else: - return None - - -def get_obs_summary_depot(t_min: Time, t_max: Time) -> MNS | None: - """ - Get observation summary from depot - """ - - jds = np.arange(t_min.jd, t_max.jd + 1) - - res = get_coverage_depot(jds) - - if len(res) == 0: - return None - - res["date"] = Time(res["obsjd"].to_numpy(), format="jd").isot - - mns = MNS(df=res) - - mns.data.query(f"obsjd >= {t_min.jd} and obsjd <= {t_max.jd}", inplace=True) - - mns.data.reset_index(inplace=True) - mns.data.drop(columns=["index"], inplace=True) - - return mns - - -def get_obs_summary(t_min, t_max=None, max_days: float = None) -> MNS | None: - """ - Get observation summary from IPAC depot - """ - now = Time.now() - - if t_max and max_days: - raise ValueError("Choose either t_max or max_days, not both") - - if t_max is None: - if max_days is None: - t_max = now - else: - t_max = t_min + (max_days * u.day) - - if t_max > now: - t_max = now - - logger.info("Getting observation logs from IPAC depot.") - mns = get_obs_summary_depot(t_min=t_min, t_max=t_max) - - if mns is not None: - logger.debug( - f"Found {len(set(mns.data['exposure_id']))} observations in total." - ) - else: - logger.debug("Found no observations on IPAC depot.") - - return mns diff --git a/nuztf/parse_nu_gcn.py b/nuztf/parse_nu_gcn.py index 7fa0262..82a8d03 100644 --- a/nuztf/parse_nu_gcn.py +++ b/nuztf/parse_nu_gcn.py @@ -180,6 +180,7 @@ def parse_gcn_circular(gcn_number: int): endpoint = f"https://gcn.nasa.gov/circulars/{gcn_number}.json" res = requests.get(endpoint) + res_json = res.json() subject = res_json.get("subject") diff --git a/nuztf/plot.py b/nuztf/plot.py index 079281f..2e2dd78 100644 --- a/nuztf/plot.py +++ b/nuztf/plot.py @@ -56,6 +56,7 @@ def lightcurve_from_alert( title: str = None, include_ulims: bool = True, include_cutouts: bool = True, + include_ps1: bool = True, include_crossmatch: bool = True, mag_range: list = None, z: float = None, @@ -124,23 +125,24 @@ def lightcurve_from_alert( ): create_stamp_plot(alert=cutout_, ax=ax_, cutout_type=type_) - img_cache = CUTOUT_CACHE_DIR.joinpath(f"{name}_PS1.png") + if include_ps1: + img_cache = CUTOUT_CACHE_DIR.joinpath(f"{name}_PS1.png") - if not img_cache.is_file(): - img = get_ps_stamp( - candidate["ra"], candidate["dec"], size=240, color=["y", "g", "i"] - ) - img.save(img_cache) + if not img_cache.is_file(): + img = get_ps_stamp( + candidate["ra"], candidate["dec"], size=240, color=["y", "g", "i"] + ) + img.save(img_cache) - else: - from PIL import Image + else: + from PIL import Image - img = Image.open(img_cache) + img = Image.open(img_cache) - cutoutps1.imshow(np.asarray(img)) - cutoutps1.set_title("PS1", fontdict={"fontsize": "small"}) - cutoutps1.set_xticks([]) - cutoutps1.set_yticks([]) + cutoutps1.imshow(np.asarray(img)) + cutoutps1.set_title("PS1", fontdict={"fontsize": "small"}) + cutoutps1.set_xticks([]) + cutoutps1.set_yticks([]) # If redshift is given, calculate absolute magnitude via luminosity distance # and plot as right axis diff --git a/nuztf/skymap_scanner.py b/nuztf/skymap_scanner.py index 44e5ba6..4d6a6fb 100644 --- a/nuztf/skymap_scanner.py +++ b/nuztf/skymap_scanner.py @@ -42,7 +42,6 @@ def __init__( self.prob_threshold = prob_threshold self.n_days = n_days self.event = event - self.rev = rev self.prob_threshold = prob_threshold self.output_nside = output_nside @@ -55,10 +54,11 @@ def __init__( self.skymap = Skymap( event=self.event, - rev=self.rev, + rev=rev, prob_threshold=self.prob_threshold, output_nside=self.output_nside, ) + self.rev = self.skymap.rev self.t_min = Time(self.skymap.t_obs, format="isot", scale="utc") @@ -428,9 +428,11 @@ def unpack_skymap(self, output_nside: None | int = None): map_coords.append(self.extract_ra_dec(nside, i)) pixel_nos.append(i) - pixel_area = hp.nside2pixarea(nside, degrees=True) * float(len(map_coords)) + total_pixel_area = hp.nside2pixarea(nside, degrees=True) * float( + len(map_coords) + ) - self.logger.info(f"Total pixel area: {pixel_area} degrees") + self.logger.info(f"Total pixel area: {total_pixel_area} degrees") map_coords = np.array( map_coords, dtype=np.dtype([("ra", float), ("dec", float)]) @@ -442,7 +444,7 @@ def unpack_skymap(self, output_nside: None | int = None): nside, self.skymap.data[self.skymap.key][mask], self.skymap.data, - pixel_area, + total_pixel_area, self.skymap.key, ) diff --git a/tests/test_observations.py b/tests/test_observations.py index 9144dbf..98a1d0f 100644 --- a/tests/test_observations.py +++ b/tests/test_observations.py @@ -7,7 +7,7 @@ from astropy import units as u from astropy.time import Time -from nuztf.observations import get_obs_summary_irsa, get_obs_summary_skyvision +from nuztf.observations import get_obs_summary_skyvision class TestCoverage(unittest.TestCase):