Merge 74d9298 into 7f39d00

.github/workflows/continous_integration.yml

@@ -58,6 +58,7 @@ jobs:
           SKYVISION_USER: ${{ secrets.skyvision_user }}
           SKYVISION_PASSWORD: ${{ secrets.skyvision_password }}
           FRITZ_TOKEN: ${{ secrets.fritz_token }}
+          DESY_CLOUD_TOKEN: ${{ secrets.desy_cloud_token }}
           ZTFDATA: ./
         run: |
           poetry run coverage run -m pytest -v

nuztf/ampel_api.py

@@ -9,12 +9,15 @@
 import backoff
 import numpy as np
 import requests
+from ampel.util.json import load
+from ampel.ztf.util.ZTFIdMapper import ZTFIdMapper
 from astropy.io import fits  # type: ignore
 from astropy.time import Time  # type: ignore
-from nuztf.credentials import load_credentials
-from nuztf.utils import deres
 from requests.auth import HTTPBasicAuth
 
+from nuztf import utils
+from nuztf.credentials import load_credentials
+
 API_BASEURL = "https://ampel.zeuthen.desy.de"
 API_ZTF_ARCHIVE_URL = API_BASEURL + "/api/ztf/archive/v3"
 API_CATALOGMATCH_URL = API_BASEURL + "/api/catalogmatch"
@@ -489,7 +492,7 @@ def ampel_api_skymap(
         gt = "$gt"
 
     # Now we reduce the query size
-    regions = deres(nside=nside, ipix=pixels)
+    regions = utils.deres(nside=nside, ipix=pixels)
 
     query = {
         "regions": regions,
@@ -663,3 +666,43 @@ def ampel_api_catalog(
         raise requests.exceptions.RequestException
 
     return res
+
+
+def get_preprocessed_results(file_basename: str, logger=None) -> None | list:
+    """
+    Access the DESY Cloud to look if there are precomputed results from an AMPEL run there
+    """
+    if logger is None:
+        logger = logging.getLogger(__name__)
+
+    desy_cloud_token = load_credentials("desy_cloud_token", token_based=True)
+
+    filename = file_basename + ".json.gz"
+
+    res = requests.get(
+        f"https://syncandshare.desy.de/public.php/webdav/{filename}",
+        headers={"X-Requested-With": "XMLHttpRequest"},
+        auth=(desy_cloud_token, "bla"),
+    )
+
+    if res.status_code != 200:
+        logger.warning(
+            "\n\n-------------------- !! -------------------\nSomething went wrong with your query.\nCheck your credentials and make sure Ampel\nhas run correctly at Desy.\n-------------------- !! -------------------\n\n"
+        )
+        return None
+
+    with open(f"{filename}", "wb") as f:
+        f.write(res.content)
+
+    res = []
+    with gzip.open(filename, "rb") as f_in:
+        data = load(f_in)
+        for t in data:
+            ztf_id = ZTFIdMapper.to_ext_id(t.stock.get("stock"))
+            pp = t.get_photopoints()
+            pp_reformatted = utils.reformat_downloaded_results(
+                photopoints=pp, ztf_id=ztf_id
+            )
+            res.append(pp_reformatted)
+
+    return res

nuztf/base_scanner.py

@@ -23,6 +23,9 @@
 from astropy.time import Time
 from gwemopt.ztf_tiling import get_quadrant_ipix
 from matplotlib.backends.backend_pdf import PdfPages
+from tqdm import tqdm
+from ztfquery import fields as ztfquery_fields
+
 from nuztf.ampel_api import (
     ampel_api_acknowledge_chunk,
     ampel_api_cone,
@@ -38,8 +41,6 @@
 from nuztf.observations import get_obs_summary
 from nuztf.plot import lightcurve_from_alert
 from nuztf.utils import cosmo
-from tqdm import tqdm
-from ztfquery import fields as ztfquery_fields
 
 DEBUG = False
 RATELIMIT_CALLS = 10
@@ -302,9 +303,10 @@ def scan_area(
         Retrieve alerts for the healpix map from AMPEL API,
         filter the candidates and create a summary
         """
-
         query_res = self.query_ampel(t_min=t_min, t_max=t_max)
 
+        ztf_ids_zero_stage = [res["objectId"] for res in query_res]
+
         ztf_ids_first_stage = []
         for res in tqdm(query_res):
             if self.filter_f_no_prv(res):
@@ -504,6 +506,7 @@ def extract_npix(nside, ra, dec):
     def create_candidate_summary(self, outfile=None):
         """Create pdf with lightcurve plots of all candidates"""
         if len(self.cache.items()) == 0:
+            self.logger.info("No candidates found, skipping pdf creation")
             return
 
         if outfile is None:
@@ -530,6 +533,10 @@ def create_candidate_summary(self, outfile=None):
 
     def create_overview_table(self, outfile=None):
         """Create csv table of all candidates"""
+        if len(self.cache.items()) == 0:
+            self.logger.info("No candidates found, skipping csv creation")
+            return
+
         if outfile is None:
             csv_path = self.summary_path + ".csv"
         else:

nuztf/credentials.py

@@ -61,7 +61,18 @@ def load_credentials(name: str, token_based: bool = False):
         logging.info('Set up "tns_api_token" credentials')
 
 except KeyError:
-    logging.info("No Token for TNS API found in environment" "Assume they are set.")
+    logging.info("No Token for TNS API found in environment" "Assume it is set.")
+
+try:
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=UserWarning)
+        io.set_account(
+            "desy_cloud_token", token=os.environ["DESY_CLOUD_TOKEN"], token_based=True
+        )
+        logging.info('Set up "desy_cloud_token" credentials')
+
+except KeyError:
+    logging.info("No Token for DESY Cloud found in environment" "Assume it is set.")
 
 try:
     with warnings.catch_warnings():
@@ -70,4 +81,4 @@ def load_credentials(name: str, token_based: bool = False):
         logging.info('Set up "fritz" credentials')
 
 except KeyError:
-    logging.info("No Token for Fritz API found in environment" "Assume they are set.")
+    logging.info("No Token for Fritz API found in environment" "Assume it is set.")

nuztf/skymap.py

@@ -207,6 +207,8 @@ def get_gw_skymap(self, event_name: str, rev: int):
         elif rev > len(voevents):
             raise Exception("Revision {0} not found".format(rev))
 
+        self.rev = rev
+
         latest_voevent = voevents[rev - 1]
         self.logger.info(f"Found voevent {latest_voevent['filename']}")
 

nuztf/skymap_scanner.py

@@ -12,15 +12,17 @@
 import yaml
 from astropy.time import Time
 from astropy_healpix import HEALPix
+from tqdm import tqdm
+from ztfquery.io import LOCALSOURCE
+
 from nuztf.ampel_api import (
     ampel_api_acknowledge_chunk,
     ampel_api_lightcurve,
     ampel_api_skymap,
+    get_preprocessed_results,
 )
 from nuztf.base_scanner import BaseScanner
 from nuztf.skymap import Skymap
-from tqdm import tqdm
-from ztfquery.io import LOCALSOURCE
 
 
 class RetractionError(Exception):
@@ -84,6 +86,30 @@ def get_full_name(self):
         else:
             return "?????"
 
+    def download_results(self):
+        """
+        Retrieve computed results from the DESY cloud
+        """
+        self.logger.info("Retrieving results from the DESY cloud")
+        file_basename = f"{self.skymap.event_name}_{self.skymap.rev}"
+
+        res = get_preprocessed_results(file_basename=file_basename)
+
+        if res is None:
+            final_objects = []
+        else:
+            final_objects = [alert["objectId"] for alert in res]
+            for alert in res:
+                self.cache[alert["objectId"]] = alert
+
+        final_objects = self.remove_duplicates(final_objects)
+
+        self.logger.info(
+            f"Retrieved {len(final_objects)} final objects for event {self.skymap.event_name} / map revision {self.skymap.rev} from DESY cloud."
+        )
+
+        self.final_candidates = final_objects
+
     def get_alerts(self):
         """Scan the skymap area and get ZTF transients"""
         self.logger.info("Commencing skymap scan")
@@ -496,25 +522,14 @@ def plot_coverage(self, plot_candidates: bool = True):
 
         return fig, message
 
-    # def interpolate_map(self, ra_deg, dec_deg):
-    #     """ """
-    #     interpol_map = self.skymap.hpm.interpolate_bilinear_skycoord(
-    #         SkyCoord(ra_deg * u.deg, dec_deg * u.deg), self.data[self.key]
-    #     )
-    #     return interpol_map
-
-    # def in_contour(self, ra_deg, dec_deg):
-    #     """ """
-    #     return self.interpolate_map(ra_deg, dec_deg) > self.skymap.pixel_threshold
-
 
 if __name__ == "__main__":
     import logging
 
     logger = logging.getLogger(__name__)
     logger.setLevel(logging.INFO)
 
-    scanner = SkymapScanner(logger=logger)
+    scanner = SkymapScanner()
     scanner.plot_skymap()
     scanner.get_alerts()
     scanner.filter_alerts()

nuztf/utils.py

@@ -8,6 +8,7 @@
 from collections import OrderedDict, defaultdict
 from json import JSONDecodeError
 
+import numpy as np
 import requests
 from astropy.cosmology import FlatLambdaCDM
 from nuztf.credentials import load_credentials
@@ -69,6 +70,33 @@ def is_tns_name(name: str) -> bool:
     return matches
 
 
+def reformat_downloaded_results(photopoints: list, ztf_id: str) -> dict:
+    """
+    Massage the TransientView photopoint output so it matches what the archive DB returns
+    """
+
+    # Find the index of the latest detection
+    latest_jd_index = np.argmax(np.asarray([pp["body"]["jd"] for pp in photopoints]))
+
+    resdict = {
+        "candid": photopoints[latest_jd_index]["id"],
+        "objectId": ztf_id,
+        "schemavsn": "3.3",
+        "publisher": "Ampel",
+        "candidate": photopoints[latest_jd_index]["body"],
+    }
+
+    prv_candidates = []
+
+    for i, pp in enumerate(photopoints):
+        if i != latest_jd_index:  # dict is readonly, we can't pop
+            prv_candidates.append(pp["body"])
+
+    resdict["prv_candidates"] = prv_candidates
+
+    return resdict
+
+
 def deres(nside, ipix, min_nside=1):
     """
     Originally from Ampel-ZTF-archive/ampel/ztf/archive/server (by JvS)

tests/test_skymap_scanner.py

@@ -115,3 +115,25 @@ def test_grb_scan(self):
         true_gcn = "Astronomer Name (Institute of Somewhere), ............. report,\n\nOn behalf of the Zwicky Transient Facility (ZTF) and Global Relay of Observatories Watching Transients Happen (GROWTH) collaborations: \n\nAs part of the ZTF neutrino follow up program (Stein et al. 2022), we observed the localization region of the GRB210927A with the Palomar 48-inch telescope, equipped with the 47 square degree ZTF camera (Bellm et al. 2019, Graham et al. 2019). We started observations in the g- and r-band beginning at 2021-09-27 02:44 UTC, approximately 2.1 hours after event time. We covered 0.1% (1.0 sq deg) of the reported localization region. This estimate accounts for chip gaps. Each exposure was 30s with a typical depth of 20.5 mag. \n \nThe images were processed in real-time through the ZTF reduction and image subtraction pipelines at IPAC to search for potential counterparts (Masci et al. 2019). AMPEL (Nordin et al. 2019, Stein et al. 2021) was used to search the alerts database for candidates. We reject stellar sources (Tachibana and Miller 2018) and moving objects, and apply machine learning algorithms (Mahabal et al. 2019) , and removing candidates with history of variability prior to the merger time. We are left with the following high-significance transient candidates by our pipeline, all lying within the 90.0% localization of the skymap.\n\n+--------------------------------------------------------------------------------+\n| ZTF Name     | IAU Name  | RA (deg)    | DEC (deg)   | Filter | Mag   | MagErr |\n+--------------------------------------------------------------------------------+\n| ZTF21acdvtxc |  -------  | 250.2336698 | +05.3908972 | g      | 21.80 | 0.21   | (MORE THAN ONE DAY SINCE SECOND DETECTION) \n| ZTF21acdvtxp |  -------  | 250.4636648 | +01.8436867 | g      | 21.33 | 0.18   | (MORE THAN ONE DAY SINCE SECOND DETECTION) \n| ZTF21acdvuzf |  -------  | 241.8979602 | +19.0755373 | g      | 20.82 | 0.17   | (MORE THAN ONE DAY SINCE SECOND DETECTION) \n+--------------------------------------------------------------------------------+\n\n \n\nAmongst our candidates, \n\nZTF21acdvtxc had upper limit problems. PLEASE FILL IN NUMBERS BY HAND!!! WISEA J164056.10+052327.1 ['UvS'-type source (0.03 arsec)]\nZTF21acdvtxp had upper limit problems. PLEASE FILL IN NUMBERS BY HAND!!! [MILLIQUAS: SDSS J164151.27+015037.0 - Likely QSO (prob = 95.0%) (0.21 arsec)]\nZTF21acdvuzf had upper limit problems. PLEASE FILL IN NUMBERS BY HAND!!! \n\n\nZTF and GROWTH are worldwide collaborations comprising Caltech, USA; IPAC, USA; WIS, Israel; OKC, Sweden; JSI/UMd, USA; DESY, Germany; TANGO, Taiwan; UW Milwaukee, USA; LANL, USA; TCD, Ireland; IN2P3, France.\n\nGROWTH acknowledges generous support of the NSF under PIRE Grant No 1545949.\nAlert distribution service provided by DIRAC@UW (Patterson et al. 2019).\nAlert database searches are done by AMPEL (Nordin et al. 2019).\nAlert filtering is performed with the nuztf (Stein et al. 2021, https://github.com/desy-multimessenger/nuztf).\n"
 
         self.assertEqual(res, true_gcn)
+
+    def test_gw_scan_desy_download(self):
+        self.logger.info("\n\n Testing GW result download from DESY \n\n")
+
+        gw_name = "S200115j"
+        prob_threshold = 0.9
+
+        self.logger.info(f"Scanning with GW {gw_name}")
+
+        scanner = SkymapScanner(
+            event=gw_name,
+            prob_threshold=prob_threshold,
+            n_days=2,
+        )
+
+        scanner.download_results()
+
+        n_expected_candidates = 122
+
+        n_retrieved_candidates = len(scanner.final_candidates)
+
+        self.assertEqual(n_retrieved_candidates, n_expected_candidates)